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I.  INTRODUCTION 


A.  BACKGROUND 

The  NASA  Orion  Crew  Exploration  Vehicle  (CEV)  program  has  more  stringent 
safety  requirements  as  compared  to  the  Apollo  missions  of  the  1960s  and  1970s.  One 
such  requirement  is  for  “anytime  return”  of  the  spacecraft  in  the  event  of  a  mission  abort 
scenario.  The  increased  safety  margins  that  need  to  be  incorporated  into  the  spacecraft 
design  begin  to  form  a  competing  set  of  requirements.  The  design  of  spacecraft,  whether 
manned  or  unmanned,  is  primarily  driven  by  the  mass  of  fuel  carried  onboard  to  perform 
maneuvers  while  performing  its  mission.  This,  in  turn,  drives  the  total  mass  of  the 
spacecraft,  which  further  constrains  the  design  for  mission  essential  equipment  and  can 
heavily  influence  the  design  and/or  selection  of  the  launch  vehicle.  Excess  safety 
margins  should  be  avoided,  however,  in  order  to  maximize  mission  effectiveness. 

The  application  of  optimal  control  theory  to  astrodynamics  applications  allows  for 
designers  to  generate  engineering  solutions  that  maximize  desirable  perfonnance 
characteristics  as  well  as  minimizing  the  need  for  excess  design  margins.  Current  design 
practices  may  unwittingly  limit  the  solutions  being  created  by  enforcing  seemingly 
practical  constraints.  As  the  spacecraft  design  matures,  clearly  there  is  a  need  to  start 
imposing  design  constraints  due  to  the  available  technology,  product  specifications,  and 
material.  However,  early  in  the  design  it  is  possible  to  find  “best”  solutions  by  limiting 
the  number  of  constraints  imposed  on  the  design  and  mathematically  validating  the 
solutions  using  existing  optimal  control  theory.  The  tools  for  solving  complex  optimal 
control  problems  allow  for  rapid  and  robust  solutions  and  straightforward  methods  are 
available  for  ascertaining  their  feasibility. 

B.  SCOPE  OF  THE  RESEARCH 

This  thesis  investigated  the  feasibility  of  using  solely  the  lower  thrust  auxiliary 
engines  of  the  Orion  CEV  for  return  to  earth  from  low  lunar  orbit.  Using  DIDO©,  which 
utilizes  MATLAB  to  solve  complex,  non-linear  optimal  control  problems,  solutions  were 
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generated  which  demonstrate  the  existence  of  trajectories  which  would  allow  the 
spacecraft  to  return  using  only  the  auxiliary  engines  while  meeting  specific  mission 
timeline  and  fuel  consumption  requirements.  Furthermore,  the  solutions  generated  are 
shown  to  meet  the  necessary  conditions  of  optimality  by  applying  Pontryagin’s  Theory. 
While  there  exists  a  range  of  possible  solutions,  two  significant  solutions  are  presented: 

(1)  the  absolute  minimum  time  required  for  the  maneuver  with  no  constraint  on  fuel,  and 

(2)  the  minimum  fuel  required  to  reach  the  given  final  boundary  conditions  within  48 
hours. 

The  first  case  is  shown  primarily  to  show  the  absolute  lower  limit  on  time 
required  for  the  maneuver  with  the  given  engine  characteristics  and  the  dynamics  of  the 
problem.  The  necessary  conditions  for  optimality  of  this  solution  are  shown  and  a 
feasibility  analysis  is  conducted  to  check  the  solutions.  The  feasibility  of  the  control 
solution  is  conducted  by  the  application  of  the  Bellman  technique.  While  the  Bellman 
technique  is  not  shown  for  the  second,  minimum  fuel  case,  it  demonstrates  a  method  to 
verify  the  fidelity  of  the  generated  solutions  and  determine  if  the  solution  is  “flyable” 
given  a  nominal  set  of  discrete  control  points.  The  second  case  shows  that  the  auxiliary 
engines  can  be  used  to  complete  the  entire  moon-to-earth  Trans-earth  Injection  (TEI) 
maneuver  sequence  and  still  meet  the  48-hour  time  requirement  and  use  less  fuel  that  the 
capacity  of  the  Orion’s  fuel  tank.  In  all  cases,  the  necessary  conditions  for  optimality  are 
demonstrated  using  Pontryagin’s  Hamiltonian  Minimization  Condition  (HMC)  and  where 
singular  arcs  appear  in  the  trajectory,  further  analysis  is  performed  to  show  optimality. 

Another  problem  investigated  in  this  thesis  is  the  control  of  the  singular  arc 
segment  of  a  previously  determined  DIDO©  solution  using  the  CEV  main  engines.  The 
main  engine  solution  resulted  in  the  return  trajectory  consisting  of  three  separate  TEI 
bums.  The  first  and  third  burns  are  maximum  thrust  burns,  known  as  bang-bang 
maneuvers.  The  first  burn  raises  the  apoapsis  of  the  low  lunar  orbit,  while  the  final  bum 
provides  the  necessary  velocity  to  escape  the  moon’s  gravity  on  its  way  towards  earth. 
The  middle  bum  is  given  as  a  finite  burn  of  varying  thrust  where  the  spacecraft  conducts 
a  plane  change  maneuver.  A  study  of  the  high  thrust  singular  arc  solution  was  conducted 
in  order  to  attain  more  in-depth  details  of  the  control  trajectory  in  order  to  determine  the 
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feasibility  of  using  the  auxiliary  engines  to  conduct  the  singular  burn  segment  of  the  high 
thrust  return  trajectory.  Figure  1  shows  the  thrust  profile  for  the  minimum  fuel,  fixed 
time  return  solution  using  the  Orion  CEV  main  engines  found  by  Yan  et  al.  (Yan,  Gong, 
Park,  Ross,  &  Souza,  2010). 
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II.  EQUATIONS  OF  MOTION 


A.  INTRODUCTION 

For  interplanetary  trajectories,  the  two  body  problem  no  longer  suffices  as  the 
spacecraft  will  pass  through  areas  of  varying  gravitational  influence  as  it  gets  further 
away  from  one  body  and  gets  closer  to  another.  Even  outside  the  sphere  of  influence  of  a 
particular  body,  the  gravitational  acceleration  of  that  body  may  provide  a  minimal,  but 
influential  gravitational  affect,  often  referred  to  as  a  gravitational  perturbation.  This 
chapter  first  defines  the  specific  reference  frame  used  for  each  case  in  this  thesis. 
Secondly,  it  describes  the  gathering  of  ephemerides  data  on  the  gravitational  bodies 
which  have  significant  effect  on  the  spacecraft  trajectory  which  is  ultimately  used  in 
setting  up  the  dynamics  of  the  moon-earth  return  mission.  Finally,  a  brief  description  of 
the  dynamics  is  given,  followed  by  an  example  of  why  the  problems  require  n=4  for 
moon-earth  return  missions. 

B.  DEFINING  THE  REFERENCE  FRAME  (J2000) 

If  one  has  to  move  from  one  specific  location  to  another,  one  has  to  have  good 
knowledge  of  the  current  position  as  well  as  an  accurate  description  or  knowledge  of  the 
intended  final  position.  In  common  tenns,  it  is  called  navigation.  On  the  surface  of  the 
earth,  the  art  of  navigation  is  generally  simple.  Barring  any  obstacles,  one  simply  takes 
the  shortest  distance  between  two  points  (which  are  fixed)  and  travels  along  the  line  that 
connects  the  two  points.  This  is  done  by  establishing  a  reference  frame.  Assuming  that 
the  travel  is  in  only  two  dimensions,  a  planar  reference  frame  can  be  used.  Longer 
distances  on  the  surface  of  the  earth  can  be  referenced  to  the  surface  of  a  sphere.  Despite 
the  rotations  of  the  earth,  the  two  dimensional  traveler  need  only  stay  on  the  line 
connecting  the  two  points.  The  reference  frame,  whether  a  flat  plane  or  surface  of  the 
sphere  is  considered  inertial  since  it  is  not  accelerating  with  respect  to  the  traveler. 

Depending  on  the  specific  application,  such  as  interplanetary  space  missions  or 
tracking  satellites  in  orbit,  a  number  of  different  reference  frames  are  used.  For  rough 
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calculations,  the  reference  systems  used  can  use  an  object’s  center  of  mass  as  the  center. 
However,  for  spaceflight  in  the  vicinity  of  the  earth  and  moon,  it  might  be  useful  to  use 
the  earth-moon  barycenter,  or  system  center  of  mass.  Synodic  systems  are  also  used 
which  are  systems  that  rotate  about  the  barycenter.  In  synodic  coordinate  systems,  it  is 
assumed  that  the  rotation  about  the  barycenter  is  with  constant  velocity  with  respect  to  an 
inertial  frame. 

The  Heliocentric  Coordinate  System  (XYZ)  is  a  sun  centered  system  with  the 
primary  axis,  X,  pointed  at  the  vernal  equinox.  It  is  a  right  handed  system  with  respect  to 
the  direction  of  earth’s  orbit  around  the  sun  with  Z  perpendicular  to  the  ecliptic  plane. 
With  the  correct  transformations,  this  coordinate  system  could  be  converted  into  a 
synodic  coordinate  system,  by  using  the  barycenter  of  the  solar  system  which  is  a  point 
slightly  off  from  the  geometric  center  of  the  sun. 


Figure  2.  Heliocentric  Reference  Frame  (From  Vallado,  1997) 


A  commonly  used  coordinate  system  in  astrodynamics  is  called  the  Geocentric 
Equatorial  Coordinate  System  (IJK)  which  is  a  non-rotating  system  with  the  primary 
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axis,  /,  pointed  at  the  vernal  equinox.  The  J  axis  is  90  degrees  to  the  east  and  the  K  axis 
points  out  of  the  earth’s  north  pole.  This  system  is  also  known  as  the  Earth  Centered 
Inertial  (ECI)  system.  The  ECI  system  in  actuality  does  move  over  time  due  to  effects 
such  as  precession  and  nutation.  The  J2000  system  is  essentially  the  ECI  system  at  a 
fixed  epoch,  previously  noted.  This  is  considered  an  inertial  frame  from  which 
calculations  known  as  reduction  formulas  can  be  used  to  advance  (or  regress)  to  another 
epoch. 


Figure  3.  Earth  Centered  Inertial  Reference  Frame  (From  Vallado,  1997) 

The  problem  with  space  travel  is  in  determining  and  characterizing  an  inertial 
reference  frame,  since  by  nature  all  the  objects  in  the  solar  system  are  in  a  state  of 
continual  motion,  and  in  most  cases  non-uniform  motion  due  to  the  effects  of 
gravitational  forces  between  them.  For  short  time  spans,  it  may  suffice  to  assume  an 
inertial  frame  with  the  earth  at  its  center;  however  the  nonuniform  motion  of  the  earth 
and  other  solar  system  bodies  have  an  increasingly  perturbing  effect  over  time. 
Therefore,  to  accurately  determine  the  position  of  a  body  (celestial  body  or  spacecraft) 
traveling  through  the  solar  system,  one  has  to  have  a  method  to  accurately  determine  the 
locations  of  all  the  gravitational  bodies,  or  at  least  the  ones  which  have  the  most 
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significant  effect  on  the  motion  of  that  body.  The  positions  of  the  solar  system  bodies 
can  be  detennined  by  observation  and  the  future  positions  can  be  determined,  assuming 
the  motions  are  well  characterized,  at  a  future  time. 

For  consistency,  it  has  been  useful  to  designate  a  specific  time,  or  epoch,  as  a 
reference  time  with  which  the  reference  frame  used  for  navigation  purposes  can  be 
advanced  or  rotated.  The  current  system  in  use  today  is  called  J2000.  The  J2000  epoch 
is  specified  to  be,  January  1,  2000  at  11:59:27.816  TAI  (International  Atomic  Time),  or 
about  noon  UT  (Universal  Time).  Therefore,  J2000  refers  to  a  specific  epoch  from  which 
to  attach  a  coordinate  system  as  a  basis  for  all  celestial  observations  (Vallado,  1997). 
Typically,  an  ECI  reference  frame  is  used  as  the  basis.  If  a  moon-centered  inertial  frame 
was  desired,  then  with  the  proper  transfonnation,  the  frame  could  be  translated  from  the 
earth  to  the  moon’s  location  at  the  J2000  epoch. 

In  addition  to  planetary  motion  about  the  sun,  the  Earth’s  position  and  orientation 
with  respect  to  the  J2000  ECI  frame  varies  due  to  gravitational  interactions  with  other 
bodies  in  the  solar  system,  but  primarily  due  to  the  Moon  and  Sun.  The  four  primary 
transformations  which  must  be  made  are  for  precession,  nutation,  sidereal  time,  and 
polar  motion.  Precession  results  from  perturbations  from  the  Sun,  Moon,  and  planets  and 
it  changes  the  orientation  of  the  ecliptic  with  respect  to  the  J2000  frame.  The  reduction 
formulas  for  precession  allow  for  the  calculation  of  the  mean  equator  of  date  from  the 
mean  equator  at  epoch  (J2000).  Nutation  is  largely  caused  by  the  moon  and  is  the 
lengthiest  transfonnation,  consisting  of  over  100  trigonometric  tenns  and  is  a  periodic 
perturbation.  Taking  into  account  the  effects  of  nutation  transforms  the  mean  equator  of 
date  into  the  true  equator  of  date.  The  third  transformation,  sidereal  time,  transforms  the 
non-rotating  true  of  date  frame  to  the  Earth-fixed  coordinate  system.  Finally,  polar 
motion  accounts  for  the  changing  location  of  the  North  Pole.  The  motion  of  the  North 
Pole  follows  a  circular  spiral  pattern  and  has  a  maximum  variation  of  about  only  9  meters 
in  any  direction  (Vallado,  1997). 
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Figure  4.  Precession  and  Nutation  (From  Vallado,  1997) 

The  problems  studied  in  this  thesis  use  a  J2000  moon-centered  frame  of  reference 
since  the  spacecraft  originates  in  close  proximity  to  the  moon.  In  this  way,  the  moon  is 
the  primary  gravitational  body  exerting  a  force  on  the  spacecraft,  while  other  celestial 
bodies,  including  the  earth,  exert  a  perturbing  force  on  the  spacecraft. 

C.  JPL  HORIZONS  DATA 

The  Jet  Propulsion  Lab  (JPL)  HORIZONS  on-line  system  provides  accurate 
ephemeris  data  for  solar  system  objects  to  include  538186  asteroids,  3066  comets,  170 
planetary  satellites,  8  planets,  the  Sun,  LI,  L2,  select  spacecraft,  and  system  barycenters 
(Jet  Propulsion  Laboratory).  Using  the  web  interface,  ephemeris  data  was  collected  for 
earth,  sun,  Mars,  and  Venus  as  target  bodies  with  the  origin  body  as  the  moon.  The 
ephemeris  data  was  collected  in  the  form  of  position  (X,  Y,  Z)  and  velocity  (VX,  VY, 
VZ)  vectors  from  the  moon  to  the  target  bodies  in  the  J2000  frame  covering  a  time  span 
of  ten  days  starting  at  midnight  on  April  2,  2024,  until  midnight  on  April  12,  2024,  in 
increments  of  one  minute.  The  vectors  were  with  respect  to  the  ecliptic  and  mean 
equinox  of  the  J2000  epoch. 
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In  order  to  have  access  to  accurate  position  vectors  to  the  target  bodies  at  any 
time,  the  relative  motion  of  the  target  bodies  can  be  characterized  as  a  smooth  function  of 
time.  In  this  study  the  Matlab  function,  polyfit,  was  used  to  be  able  to  compare  the 
accuracy  of  the  data  used  by  Yan  et  ah,  who  used  the  same  function  (Yan  et  ah,  2010). 
Even  though  the  data  points  are  very  close  together  in  comparison  to  the  total  duration,  it 
is  desired  to  create  a  smooth  function  in  order  to  determine  the  position  vectors  between 
the  discrete  points.  The  goal  for  the  overall  project  is  to  ensure  that  the  accuracy  of  the 
position  of  the  spacecraft  and  the  celestial  bodies  is  within  one  kilometer.  At  spacecraft 
interplanetary  velocities,  one  minute  is  actually  a  long  duration.  Take  for  example  a 
spacecraft  travelling  at  a  velocity  of  7.7  km/sec.  In  one  minute,  the  spacecraft  has 
travelled  462  km. 

The  polyfit  function  in  Matlab  interpolates  between  the  data  points  by  fitting  a 
polynomial  of  specified  order  between  the  data  points,  x,  by  solving  for  P(x)  in  a  least 
squares  approximation.  In  order  to  minimize  calculation  errors  when  solving  for  the  state 
vectors  using  polynomials,  the  order  of  the  polynomial  should  be  minimized  but  balanced 
against  desired  accuracy.  For  the  position  vectors,  an  accuracy  on  the  order  of  one  meter 
was  desired  and  was  accomplished  by  using  polynomials  of  degree  n  =  12.  The 
coefficients  for  the  position  and  velocity  vector  polynomials  are  given  in  TABLES  1  and 
2.  The  position  and  velocities  of  the  bodies  can  then  be  determined  simply  by  applying 
the  following  equation,  where  P  applies  to  either  the  Earth  or  Sun  position  or  velocity  at 
the  specific  time. 

n 

P  =  ^f^akt"  k  ,t  =  days  from  start  of  epoch  ( 1 . 1 ) 

k= 0 

Figures  5  and  6  show  the  polynomial  fitting  errors  for  the  sun  and  earth  positions 
taken  from  the  Horizons  ephemeris  database;  increasing  the  order  of  the  polynomial  did 
not  result  in  an  increase  in  accuracy. 
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Table  2.  Moon-to-Earth  Coefficients 


12 


Error  (km) 


0.015 


Polyfit  Accuracy 


0.01 


0.005 


-0.005 


-0.01 


-0.015 


2  4  6  8 

Julian  Day  (JD  =2455347.5) 


Figure  5.  Sun  Position  Error  wrt  Moon 
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Figure  6.  Earth  Position  Error  wrt  Moon 
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D. 


THE  IMPORTANCE  OF  ACCURACY 


One  might  ask  why  it  is  so  important  to  take  into  account  what  seems  like  very 
small  variations  and  perturbations.  Obviously,  over  long  periods  of  time,  the  effects  of 
the  perturbations  would  increase,  but  these  could  be  taken  into  account  using  frequent 
observations  and  updates  to  spacecraft  navigation  systems  and  simply  make  course 
adjustments  during  flight.  The  greatest  constraint  on  spacecraft  is  typically  its  mass,  due 
to  launch  constraints  imposed  by  launch  vehicles.  A  portion  of  this  mass  is  fuel  which 
the  spacecraft  requires  for  orbit  transfers,  station  keeping,  and  attitude  control.  While  the 
spacecraft  could  conceivably  make  correction  maneuvers  in  flight,  it  comes  at  a  cost  of 
fuel.  Any  additional  fuel  to  a  spacecraft  comes  at  a  cost  to  allocated  mass  in  other 
systems.  For  interplanetary  missions,  especially  for  future  manned  exploration  missions, 
the  payload  and  support  systems  are  critical  to  maximize  mission  duration  and  ultimately 
mission  viability.  In  practice,  the  spacecraft  design  should  minimize  the  amount  of  fuel 
required  to  conduct  its  mission.  Therefore,  for  design  purposes,  it  is  important  to  take 
into  account  additional  gravitational  effects  that  act  on  the  spacecraft. 

If  the  target  body’s  motion,  such  as  the  Moon  or  Mars,  is  very  accurately 
determined,  then  optimal  trajectories  and  thrust  regimes  could  be  used  without  the  need 
for  correction  maneuvers  during  transit.  An  analogy  might  be  a  bullet  travelling 
downrange  to  a  target  at  a  distance  of  1,000  meters.  If  the  bullet’s  initial  trajectory  is  off 
by  0.1  degrees,  it  will  only  miss  the  center  of  the  target  by  1  millimeter.  However, 
suppose  a  spacecraft  travelling  in  a  straight  line,  travels  a  distance  of  80,000,000 
kilometers,  the  rough  straight  line  distance  between  Earth’s  orbit  around  the  Sun  and 
Mars’  orbit.  A  trajectory  error  of  the  same,  0.1  degrees  will  have  the  spacecraft  miss  its 
intended  target  by  over  120  kilometers.  Considering  the  great  distance  travelled,  this 
may  not  seem  like  a  significant  number,  however  when  trying  to  hit  the  target  on  the 
mark,  with  minimum  fuel,  this  is  indeed  a  considerable  error.  This  highlights  the  need 
for  a  very  accurate  description  of  the  motion  of  the  planetary  bodies  and  this  primarily 
relies  on  an  accurate  accounting  of  time  from  the  inertial  reference  frame  time. 
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For  a  spacecraft  returning  to  Earth  from  the  Moon,  in  order  to  minimize  fuel, 
optimal  control  trajectories  are  being  created  such  that  the  spacecraft  has  a  1  km  window 
for  earth  orbit  insertion.  This  means  that  in  order  for  the  trajectory  to  be  accurately 
designed,  the  spacecraft  dynamics  need  to  be  of  very  high  order  of  accuracy.  It  can  be 
shown  that  for  a  moon-earth  return  trajectory,  the  perturbation  effects  of  the  Earth  and 
Sun  must  be  taken  into  account  in  order  to  achieve  the  accuracy  required.  In  order  to 
demonstrate  this  fact,  a  portion  of  the  trajectory  from  Yan  et  al.  moon  to  earth  trajectory 
using  the  CEV  main  engines  was  investigated  to  show  the  effects  of  changing  the  number 
of  bodies  included  in  the  dynamics  equations  (Yan  et  al.,  2010).  The  segment 
investigated  has  the  initial  position  immediately  following  the  first  TEI  maneuver  and  the 
target  position  immediately  prior  to  the  second  TEI  maneuver.  The  trajectory  was 
propagated  using  the  MatLab  ode45  function  and  allowed  to  propagate  until  the  time  of 
the  second  burn  maneuver.  Starting  with  a  simple  restricted  two-body  dynamics  equation 
using  only  the  Moon’s  gravity  (Body  2),  successive  runs  were  conducted,  each  adding  the 
gravitational  effects  of  the  Earth  (Body  3),  Sun  (Body  4),  and  finally  Mars  (Body  5).  The 
time-varying  positions  of  the  earth  and  sun  with  respect  to  the  moon  were  taken  into 
account  as  described  previously  as  well  as  the  governing  equations  of  motion  as 
described  in  the  following  section.  The  initial  conditions  are  as  follows: 
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km 

To 

255.0671214283607 

km 

zo 
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Figures  7  and  8  show  that  the  four-body  dynamics,  which  includes  the  moon, 
earth  and  sun,  satisfy  the  accuracy  requirements  while  the  addition  of  the  fifth  body  do 
not  contribute  a  significant  effect  on  the  trajectory,  therefore  its  effect  need  not  be 
included. 
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Figure  7.  Effects  of  Additional  Gravitational  Bodies  on  Displacement 


Figure  8.  Effects  of  Additional  Gravitational  Bodies  on  Velocity 
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E.  EQUATIONS  OF  MOTION 


1.  The  Restricted  Two-Body  Problem 

In  the  restricted  two-body  dynamic  model,  the  mass  of  the  spacecraft  is  deemed 
negligible  and  equation  of  motion  is: 


r 


GMfr 3 


\rj 


(1.2) 


This  equation  describes  the  force  of  attraction  between  two  bodies,  namely  a 
planet  (moon  or  sun)  and  a  spacecraft  where  G  is  the  universal  gravitational  constant,  M 
is  the  mass  of  the  large  body,  m  is  the  mass  of  the  spacecraft,  and  r  is  the  position  vector 
indicated  by  the  negative  sign  to  be  from  the  spacecraft  to  the  large  body  (Figure  9).  For 
generality,  the  origin  of  the  inertial  reference  frame  will  not  be  located  at  the  center  of  the 
large  mass.  Therefore,  the  position  vector  r  will  be  the  difference  between  the  vectors 
from  the  origin  to  each  of  the  bodies  as  follows: 


r=rsc~rM 


(1.3) 
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Figure  9.  2-body  Dynamic  Model 

Since  the  data  from  HORIZONS  is  in  XYZ  vector  format,  the  r  vector  is 
converted  into  Cartesian  coordinates,  substituting  // M  for  GM  .  Three  equations  of 
motion  are  generated,  one  for  each  axis  of  motion  where: 


x  =  — 


Mmx 


2  2 

+  Z 


y- 


z  =  — 


(x2+y2+z  ^y]x  +y 

_ M _ 

(x2+y2+z2yx2+y2+z2 

_ EmI _ 

(x  +y  +z  Jyjx  +y  +z 


1  ,  X  XSC  XM 


j,  y  =  ysc-yM 


k,  z  zsc  zM 


(1.4) 
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2. 


The  n-Body  Problem 


From  the  point  of  view  of  an  inertial  observer,  the  restricted  n-body  equation  of 
motion  with  respect  to  the  primary  body  is  shown  in  Equation  (2.5)  below.  Figure  10 
shows  a  diagram  of  the  n-body  geometry  when  there  are  only  three  bodies  (n=3). 


ru  = ' 


/VI: 


1,2 


■2>y 

j= 3 


(  f 

rj,  2 

r3 

Vr  J.  2 


+  ■ 


r  A 

ly 

7  w7 


(1.5) 
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Figure  10.  Three-Body  Geometry  (After  Vallado,  1997) 

The  first  term  in  Equation  (1.5)  is  the  2-body  part  of  the  equation  of  motion, 
where  the  subscript  “1”  represents  the  primary  gravitational  body.  The  spacecraft  in  the 
equation  is  body  number  “2.”  The  summation  term  on  the  right  represents  the 
gravitational  perturbation  forces  due  to  additional  bodies  j  through  n,  starting  with  body 
number  “3.”  The  perturbation  term  itself  has  two  parts,  the  left  tenn  called  the  direct 
effect  because  it  represents  the  force  acting  directly  on  the  spacecraft  by  the  body.  The 
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right  hand  perturbation  term  is  called  the  indirect  effect  because  it  represents  the 
gravitational  effect  on  the  primary,  or  number  one,  body.  This  form  of  the  equation  of 
motion  for  the  n-body  problem  is  known  as  the  relative  form  and  is  used  because  the 
motion  of  the  spacecraft  is  calculated  relative  to  the  primary  body  (Vallado,  1997). 

The  indirect  effect  of  the  perturbation  effects  can  be  neglected  if  the  motion  of  the 
n-bodies  is  adequately  characterized.  Because  the  relative  positions  between  the 
gravitational  bodies  can  be  found  by  celestial  observation,  the  gravitational  forces  acting 
between  the  bodies  can  be  ignored.  Using  the  available  JPL  Horizons  data,  it  is  assumed 
that  the  effects  of  the  gravitational  bodies  on  each  other  are  already  taken  into  account, 
therefore  the  only  perturbation  effects  that  need  be  considered  acting  on  the  spacecraft 
are  the  direct  effects.  Therefore,  Equation  ( 1 .5)  reduces  to: 


r,^=-- 


/VI: 


-2A 


1,2 


j= 3 


v  r3J’2  J 


(1.6) 


3.  Equations  of  Motion  for  Moon  to  Earth  Return  Mission 

The  specific  equations  of  motion  used  in  this  thesis  are  driven  by  the 
establishment  of  a  J2000  moon-centered,  translating  Cartesian  frame.  The  origin  is  fixed 
to  the  moon’s  center  and  the  primary  axes  are  aligned  with  the  J2000  sun-centered 
inertial  Cartesian  frame  (SXYZ).  The  SXYZ  frame  has  the  X-Y  plane  aligned  with  the 
plane  of  Earth’s  orbit  and  the  X  direction  points  to  the  First  Point  of  Ares.  The  SXYZ  is 
a  right  handed  system,  with  the  Z  axis  perpendicular  to  the  X-Y  plane. 

Cartesian  coordinates  are  selected  first,  for  ease  of  using  JPL  Horizons  ephemeris 
data,  and  secondly  because  it  avoids  singularities  in  spherical  coordinates  that  may  arise 
due  to  trigonometric  constraints  (Yan  et  ah,  2010). 
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Figure  11.  J2000  Sun  and  Moon  Centered  Frame  (From  Yan  et  al.,  2010) 

The  position  of  the  spacecraft  is  represented  by  rp/s  =  rM/s  +  fp/M  ,  where  rMIS  is 
the  position  of  the  Moon  with  respect  to  the  center  of  the  Sun,  and  rp/M  is  the  position  of 

the  spacecraft  with  respect  to  the  Moon.  Taking  into  account  that  the  position  of  the 
spacecraft  is  with  respect  to  the  moon  in  a  sun-centered  inertial  frame,  then  the  following 
equation  holds  for  rp/s : 


(1.7) 


rp,S=rM,S+Xex+yey+Zez 


Defining  the  velocity  as: 


(1.8) 


Taking  the  first  and  second  derivative  of  Equation  (1.7)  with  respect  to  the  J2000 
Sun-centered  inertial  frame  yields: 
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(1.9) 


rpis  =  rM,s  +  vxex  +  vyey  +  vzez  +  a>ws  x  rP/M 


rp/S  rM/S  +  Vxex  +  Vyey  Vz^z  +  ® MIS  X  rP!M  ® M/S  X  rP!M 


(1.10) 


Noting  that  o)M/s  =  0  because  the  orientation  of  the  Moon-centered  frame  is 
aligned  with  the  Sun-centered  frame,  Equation  (1.10)  simplifies  to: 


rp/s 


=  rM/S+Vxex+Vyey+Vzez 


(1.11) 


The  fully  expressed  spacecraft  dynamics  also  have  to  be  expressed  in  the  Sun 
centered  inertial  frame  and  will  include  the  dynamics  as  in  Equation  (1.11)  as  well  as  the 
controls  or  forces  that  the  spacecraft  will  generate  as  part  of  the  complete  equations  of 
motion.  There  are  three  control  variables  which  are  necessary  for  describing  the  three 
dimension  Cartesian  components  of  the  thrust.  T  is  the  magnitude  of  the  thrust  and  a  and 
P  respectively  define  the  azimuth  and  elevation  angles. 

The  dynamics  can  be  succinctly  expressed  by  applying  Newton’s  second  law  in 
vector  form  as  follows: 


mrP/s=GM+GE+Gs+T  (1.12) 

Where  m  is  the  mass  of  the  spacecraft,  T  is  the  thrust  vector,  and  GM ,  GE , 
and  Gs  are  the  gravitational  forces  of  the  Moon,  Earth,  and  Sun  respectively.  The  Moon’s 
translating  acceleration,  also  a  required  portion  of  the  overall  system  dynamics  can  be 
defined  as: 
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(1.13) 


XM 

y m 


The  thrust  vector  is  defined  in  the  Moon  centered  frame  defining  a  as  the  azimuth, 
and  P  as  the  elevation  by  the  following: 

T cos  a  cos p 

T sin  a  cos/?  (1.14) 

T  sin  p 

The  time-varying  components  of  the  displacement  between  the  spacecraft  and  the 
center  of  the  Moon  are  given  by  (x,y,z).  Between  the  perturbing  bodies  of  the  sun  and 

earth,  and  the  primary  body  (moon),  they  are  ( xs, ys, zs)  and (xE,yE,zE)  respectively  so 

that  the  magnitude  of  the  displacements  of  the  spacecraft  from  the  gravitational  bodies 
can  be  defined  by: 

rE  =  '1(x~xe)2 +{y-yE)2  +{z~ze)2  O-15) 

rs  =  \l{x~xs)2 +{y-ys)2 +{z~zs)2 

With  (/JM,jUE,fJs) representing  the  gravitational  constants  for  the  Moon,  Earth, 

and  Sun  respectively  and  ve  representing  the  exhaust  velocity  of  the  spacecraft,  the  final 

form  of  the  system  dynamics  can  be  written,  which  includes  the  time-varying  mass  of  the 
rocket  in  the  following  equations  of  motion.  Equation  (2.16)  becomes  the  governing  set 
of  equations  of  motion  that  are  used  in  the  problems  solved  in  this  thesis. 
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(1.16) 
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III.  OPTIMAL  CONTROL  THEORY 


A.  THE  OPTIMAL  CONTROL  PROBLEM 

Problems  concerning  “optimization”  have  a  long  history  dating  back  to  the 
ancient  Greeks  where  through  geometry  they  were  able  to  determine  answers  to  problems 
such  as  the  shortest  distance  between  two  points  (a  line)  and  the  maximum  area  enclosed 
by  a  given  length  of  perimeter  (a  circle).  These  two  problems  are  part  of  a  class  of 
problems  known  as  optimization  theory,  where  the  former  is  a  minimization  problem  and 
the  latter  a  maximization  problem.  The  calculus  of  variations  is  a  theory  that  emerged 
during  the  eighteenth  century  that  deals  with  the  optimization  of  integrals.  John 
Bernoulli  (1667-1748)  posed  a  problem  concerning  the  minimization  of  an  integral  in  his 
famous  brachistochrome  problem  (Greek,  brachist  =  shortest,  chromos  =  time)  of  1696 
involving  a  bead  sliding  under  gravity  along  a  smooth  wire  in  which  he  asks  what  shape 
the  wire  must  be  to  reach  the  end  point  in  minimum  time.  During  the  1950s,  the  theory 
of  optimal  control  emerged  as  a  direct  consequence  of  space  exploration  efforts  by  the 
Americans  and  the  Soviets.  The  theory  would  allow  a  spacecraft’s  trajectory  to  be 
controlled  in  such  a  manner  as  to  reach  a  desired  destination  using  minimal  fuel  and  in 
minimum  time  (Pinch,  1993). 

Optimal  control  theory,  however,  is  not  limited  to  the  physical  guidance  and 
control  of  a  spacecraft.  It  has  applications  for  mission  planning  and  in  spacecraft  design. 
Consider  the  case  for  a  manned  spacecraft  mission  in  which  time  is  a  critical  design 
factor.  The  mission  duration  drives  requirements  such  as  the  fuel  required  for  navigation 
and  control  of  the  spacecraft  and  amount  of  life  support  equipment  and  supplies  needed 
to  sustain  the  crew.  Once  the  mission  objectives  are  specified  and  the  physics,  or 
dynamics,  of  the  problem  are  understood,  the  designer  will  seek  for  feasible  solutions  to 
problems  such  as  the  minimum  time  to  be  able  to  complete  the  mission  and  how  much 
fuel  to  carry  can  be  solved.  Therefore,  optimal  control  theory  can  be  used  in  the  design 
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process  to  determine  feasible  solutions  to  engineering  problems.  The  strength  of  the 
results  gained  lie  in  the  fact  that  the  theory  has  a  sound  basis  in  mathematics,  and  can 
therefore  be  proved  in  a  straightforward  manner. 

B.  DERIVATION  OF  THE  NECESSARY  CONDITIONS  USING  THE 
CALCULUS  OF  VARIATIONS 

The  problem  being  solved  is  to  detennine  the  controls (t) ,  such  that  the  path,  or 
trajectory,  from  an  initial  state  x(i0)  to  a  final  state  x[tf )  minimizes  the  some  cost 
functional,/  [x] ,  where: 


J\x\  =  J  F(x,u)dt  + 


(2.1) 


The  integral  term  represents  the  running  cost  and  represents  the  end-point 

cost.  Since  real  systems  are  being  considered,  it  is  safe  to  assume  that  any  trajectory 
ofx(t)  will  be  continuous  and  twice  differentiable  within  the  domain  of  the  state  space. 
It  is  also  allowable  to  have  the  controls  be  piecewise  continuous  and  bounded  within 
some  finite  control  space.  In  order  to  show  that  a  curve  x  =  x(t)  is  optimal,  necessary 
conditions  can  be  shown,  under  some  assumptions  of  smoothness,  using  the  calculus  of 
variations  and  are  as  follows  (Pinch,  1993): 

Hamiltonian 

H(A,x,u )  =  F(x,u)  + AT  f(x,u ) 

wU:;U)  (22) 


Euler-La grange  Equation 


(2.3) 
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Adjoint  Equation 


(2.4) 


M‘) 


dH 

dx 


Transversality  Condition 

*('')=§  (Z5) 

Where,  E[v_,x(tf^  :=  E^x^t^  +  vje^x^t^  is  called  the  Endpoint  Lagrangian. 


The  variable,  A ,  represents  the  co-state  and  while  it  has  no  physical  meaning,  it 
has  corresponding  values  that  mirror  the  state  variables.  The  adjoint  equation  represents 
the  co-state  dynamics  and  has  a  direct  relationship  to  its  corresponding  state  variable  via 
the  Hamiltonian  equation.  The  transversality  condition  looks  at  the  endpoint  cost  based 
on  the  final  state  in  order  to  find  final  conditions  of  the  co-states. 

C.  PONTRYAGIN’S  PRINCIPLE 

In  1955,  it  was  discovered  that  a  problem  existed  with  the  Euler-Langrange 
equations  in  that  it  was  impossible  to  solve  for  bounded  controls.  The  problem  was  not 
with  the  engineering  or  the  physics  of  the  controls  in  question,  rather  with  the  math. 
Pontryagin  then  invented  a  new  theory  known  as  the  Hamiltonian  Minimization 
Condition  (HMC)  which  then  replaced  the  Euler-Lagrange  equations.  The  HMC  was  able 
to  handle  bounded  controls,  which  is  typically  the  case  in  any  control  system.  In  most 
cases,  steering  or  attitude  control  is  unbounded  since  all  the  angles  in  a  circle  can  be 
described  in  positive  or  negative  multiples  of  2 n  .  In  the  case  of  an  applied  torque  or 
thrust  for  control,  the  realities  of  physics  dictate  that  there  will  be  upper  and  lower  limits 
to  such  a  control.  Therefore,  Pontryagin’ s  HMC  was  therefore  a  necessary  development 
in  the  practice  of  optimal  control  theory  (Ross,  2009). 

An  example  of  the  constrained  control  problem  can  be  described  in  a  two- 
dimensional  case  of  an  orbit  transfer  problem  where  the  state  and  controls  are  given  by 
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x  =[x,y,vx,vY,m]T  and  u  =\T,j3]T .  The  state  is  described  by  the  position  and  velocity 

in  two  dimensions,  and  mass.  The  control  variables  describe  the  magnitude  of  thrust  and 
angle  with  respect  to  the  reference  frame.  The  Hamiltonian  for  this  problem  becomes: 


H (A,x,u)  =  Axvx  +  Av  +  Avx  Kx  ( t )  + — cos  (3 

V  m  J 


+  A. 


Kr(t)  +  —smf3 

m  j 


\  ( 
+  A_ 


v 


(2.6) 


The  portion  of  the  dynamics  equations  that  do  not  depend  on  the  controls  are  given 
by  if,  (7).  The  exhaust  velocity,  ve,  is  constant  and  is  a  function  of  the  physical 

characteristics  of  the  rocket  motor.  The  dynamics  portions  of  the  equation  are  excluded 
for  brevity  because  they  do  not  contain  any  control  components. 

From  the  Euler-Lagrange  equation,  the  optimal  steering  law  is  determined  by: 

—  =  -^sin/?  +  — cosy9  =  0  (2.7) 

dj3  m  m 

tan  j8  =  ^  (2.8) 

A* 

However,  determining  the  optimal  thrust  law  proves  problematic  since  T  vanishes: 

^-  =  —  cos/?  +— sin/?-  —  =  0  (2.9) 

GT  m  m  ve 


The  Euler-Lagrange  equations  give  areas  where  the  slope  of  the  Hamiltonian  is 
zero,  thus  giving  potential  locations  where  maxima  and  minima  may  exist.  However,  if 
the  slope  is  not  zero  at  the  endpoints,  then  these  points  would  be  excluded  even  though 
they  might  also  contain  maxima  and  minima.  Example  Hamiltonian  curves  are  given  in 
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Figure  12  where  the  Euler-Lagrange  equations  would  give  an  incorrect  value  of  the 
controls  if  the  minimized  Flamiltonian  was  only  realized  where  the  slope  was  zero. 


Figure  12.  Flamiltonian  as  a  function  of  u  only  (From  Ross,  2009) 


Pontryagin’s  HMC,  therefore  replaces  the  Euler-Lagrange  equation  and 
minimizes  H  (2,x, «)  subject  to  the  upper  and  lower  bounds  of  the  control,  uf  <u<uu. 

He  also  proved  that  the  minimized  Hamiltonian  is  constant  with  respect  to  time  and  that 
for  minimum  time  problems  the  constant  is  equal  to  zero  and  for  fixed  time,  minimum  fuel 
problems  is  equal  to  - 1 . 

It  can  be  shown  that  from  the  HMC  that  a  relationship  exists  between  the 
direction  of  the  thrust  vector  and  the  direction  of  the  co-state  related  to  velocity. 
Equation  (3.10)  shows  that  the  angle  between  them  is  180  degrees. 


T 

T 


(2.10) 
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D.  THE  SWITCHING  FUNCTION 

Noting  that  the  Hamiltonian  as  a  function  of  T  is  linear,  and  is  bounded  subject 
toO<r  <Tmax,  the  right  hand  side  of  (2.9),  which  is  the  coefficient  of  T  in  H,  becomes 
the  switching  function,  S  . 


(2.11) 


Figure  13  shows  the  Hamiltonian  as  a  linear  function  of  T,  where  S  is  either  the 
positive  or  negative  slope  of  the  line.  This  result  implies  that  when  the  switching 
function  is  negative,  the  thrust  is  Tmax  and  when  positive,  the  thrust  is  equal  to  zero. 


These  two  cases  refer  to  what  is  known  as  bang-bang  control,  where  the  control  is  either 


maximum  or  minimum.  There  is  a  special  case  where  the  switching  function  is  equal  to 
zero  which  implies  that  the  thrust  lies  between  the  maximum  and  minimum  bounds  that  is 
covered  in  the  next  section.  It  can  be  shown  that  for  unconstrained  steering  controls,  the 
HMC  yields  the  same  result  as  the  Euler-Lagrange  equation,  so  Pontryagin’s  Principle  is 
the  general  case  of  the  Euler-Lagrange  condition  (Ross,  2009). 


H  as  a  function  of  T  only 


H  as  a  function  of  T  only 


(a) 


(b) 


Figure  13.  H  as  a  function  of  T  only  (From  Ross,  2009) 
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E. 


NECESSARY  CONDITIONS  FOR  SINGULAR  ARC  OPTIMALITY 


The  case  where  the  switching  function  equals  zero  occurs  when  an  intermediate, 
or  singular,  thrust  arc  appears  as  part  of  the  optimal  trajectory  where  the  thrust  will 
assume  a  value  between  the  upper  and  lower  bounds.  Singular  optimal  control  theory 
was  developed  as  a  result  of  the  examination  of  intennediate  thrust  arcs  in  central 
inverse-square  gravity  fields  and  what  resulted  was  the  Higher  Order  Maximum 
Principle.  The  Principle  follows  the  results  of  Pontryagin’s  HMC  and  then  takes  the 
derivative  of  the  switching  function  repeatedly  to  develop  the  necessary  conditions  for  an 
optimal  singular  arc.  This  requires  the  derivation  of  an  explicit  function  of  the  thrust  in 
order  to  develop  a  singular  optimal  control  law.  Since  the  switching  function  equals  zero, 
each  subsequent  derivative  must  clearly  also  equal  zero.  Taking  the  derivative  of  the 
switching  function  with  respect  to  time  shows  that  T  does  not  appear  explicitly  until  the 
fourth  derivative,  which  implies  that  the  problem  contains  a  second  order  singular  control 
(Park,  C.,  Yan,  H.  Gong,  Q.,  and  Ross,  I.M.,  2010).  Detennining  that  the  fourth  order 
derivative  is  equal  to  zero  is  especially  important  when  verifying  numerical  solutions  as 
the  results  are  typically  close,  but  never  exactly  equal  to  zero. 

The  problem  Yan  et  al.  were  solving  was  identical  to  the  problem  solved  in  this 
thesis,  the  only  difference  in  that  the  rocket  specifications  were  for  the  Orion  spacecraft 
main  engines.  Since  the  equations  of  motion  and  dynamics  of  the  problems  are  identical, 
the  derived  necessary  conditions  hold  here. 


dS  _  y;tv 
dt 


(2.12) 


d2S 
dt 2 


(2.13) 


d3S 
dt 3 


(2.14) 
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(2.15) 


d4S  _  1  d 

dtA  m^jXv  ■  Xv  dt 

The  resulting  third  and  fourth  order  derivative  equations  are  rather  long,  and  the 
control  variable  T  only  appears  explicitly  at  the  fourth  derivative.  The  explicit  forms  for 
the  first  through  fourth  derivatives  of  S  with  respect  to  time  are  found  in  Appendix  A. 
With  these  equations,  it  is  now  possible  to  verify  the  optimality  of  singular  arcs  that 
appear  as  part  of  an  optimal  earth  return  trajectory,  however  as  alluded  to  previously, 
numerically  validating  equality  to  zero  can  be  troublesome.  Therefore  applying  another 
result  using  the  calculus  of  variations  known  as  the  Generalized  Legendre-Clebsch,  or 
General  Convexity,  condition,  see  Equation  (3.16),  the  numerical  verification  is  made 
straightforward  (Bryson,  1975).  In  the  equation  below,  the  subscript  i  represents  the 
tenns  for  the  earth  and  sun.  Terms  without  a  sub-script  are  the  terms  for  the  moon. 


dT  dt4 


dH ' 
~dT , 


d  d4S 
dT  dt 4 


>  0 


9//(r-/lv)  9//(r-T)  15//(r-/lv)3 

rs  +  r7(l,.l,) 


GW„) 


(2.16) 


F.  TRAN S VERS ALITY  CONDITION 

The  derivation  of  the  transversality  conditions  are  outlined  in  Appendix  B  and  the 
following  condition  holds  for  the  minimum  fuel  problem. 


(2.17) 
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G.  SUMMARY  OF  NECESSARY  CONDITIONS  FOR  OPTIMALITY 

In  summary,  the  following  necessary  conditions  hold  for  optimal  control: 

1.  Hamiltonian  is  constant,  (-lfor  minimum  time,  0  for  fixed-time,  minimum 
fuel) 

2.  f  = 

T  2 


3.  r  = 


0,  when  S  >  0 

rmax  when S  <0 

Singular  when  S  =  0 


>•  where  S  =  -  —  -  — 
m  v. 


4.  Verify  optimality  of  singular  arcs. 


C^4  q 

a.  Numerically  verify  that  =  0  for  singular  arcs. 


b.  Check  Generalized  Legendre-Clebsch  Condition, 


d  d4S 

dr  df 


>  0 


5 .  Am  (tf)  =  -1  (for  minimum  fuel  cases  only) 


H.  GENERIC  PROBLEM  FORMULATION 

The  format  for  the  problems  solved  in  this  thesis  follow  Pontryagin’s  method  for 
solving  optimal  control  problems.  The  format  is  useful,  as  the  fonnulation  of  the 
problems  in  DIDO©  are  made  easier  as  the  boundary  conditions  and  constraints  are 
clearly  stated.  There  are  primarily  two  types  of  optimal  controls  problems  that  are 
solved,  namely  a  minimum  time  and  minimum  fuel.  Each  type  of  problem  has  a  special 
set  of  boundary  conditions  that  are  expressed  in  the  problem  formulation.  For  example, 
the  minimum  time  solution  will  have  a  free  condition  on  the  final  mass  and  final  time.  In 
the  coding  of  the  problem,  the  final  mass  will  have  to  be  fixed  at  zero,  or  close  to  zero  to 
prevent  an  unrealistic  negative  mass  solution.  In  the  minimum  fuel  problem  formulation, 
the  final  time  would  be  fixed  and  the  final  mass  would  be  free. 
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Each  of  the  problems  solved  in  this  thesis  include  a  third  physical  dimension,  z, 
start  with  the  following  format  and  the  conditions  derived  for  the  two-dimensional  case 
described  earlier  still  hold. 


< 


Generic  Problem  Formulation 
Where  x  =  [x, y, z, vx , vy , v. , m]  and  u=[T,a,/3], 

Minimize  =  tf  or  mf 

Subjectto  X  =  f{x,u,t ) 

x  (t0)  =  x°,  x{tf\  =  xf 

And  any  other  boundary  conditions  and/or  constraints 
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IV.  USING  DIDO© 


A.  INTRODUCTION 

The  use  of  computer-based  software  packages  such  as  DIDO©  allow  for  solving 
of  complex  optimal  control  problems,  however  there  are  many  challenges  associated  with 
solving  the  types  of  problems  discussed  in  this  thesis.  Primarily,  a  solid  understanding  of 
the  underlying  physics  of  the  problem  being  solved  should  be  well  understood,  especially 
during  the  verification  and  validation  process.  This  also  includes  an  understanding  of 
what  is  being  optimized  in  the  problem  and  checking  to  see  if  the  result  makes  physical 
sense  as  well  as  meeting  the  necessary  conditions  for  optimality.  Finally,  an 
understanding  of  the  limits  of  numerical  computation  is  useful  in  order  to  make  the  code 
run  more  efficiently.  The  verification  and  validation  steps  taken  for  the  problems  in  this 
thesis  are  described,  following  a  brief  description  of  how  the  problems  were  formulated 
in  DIDO©. 

B.  PROBLEM  FORMULATION  IN  DIDO© 

DIDO©  is  a  software  module  that  works  with  Matlab.  It  implements  Legendre 
pseudospectral  methods  which  are  used  to  solve  complex,  non-linear  optimal  control 
problems  (Ross  I.M.  ,  2007).  Its  primary  advantage  is  that  it  can  solve  the  problems  with 
as  few  or  many  discrete  points,  or  LGL  nodes,  as  defined  by  the  user,  which  makes  the 
computation  time  relatively  short.  However,  this  requires  at  iterative  process  of 
increasing  the  number  of  nodes,  until  the  desired  fidelity  is  achieved.  It  has  a 
straightforward  interface  that  allows  the  user  to  set  up  an  optimal  control  problem  using 
the  previously  described  Generic  Problem  Formulation.  The  problem  formulation  is 
broken  up  into  two  parts,  namely  the  cost  function,  which  specifies  what  design 
parameter  is  being  minimized,  and  the  constraints.  The  constraints  can  be  further  divided 
into  two  parts.  The  spacecraft  in  each  of  the  problems  is  constrained  to  move  subject  to 
the  equations  of  motion,  or  the  dynamics.  In  each  of  the  problems,  there  is  also  a  set  of 
boundary  conditions  describing  the  initial  and  final  states  of  the  spacecraft  which  were 
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given  parameters.  The  problem  formulation  in  DIDO©  follows  this  separation  of  the 
problem  formulation  by  requiring  separate  M-files  that  are  called  by  a  main  problem  file. 
An  additional  file  may  be  used  to  specify  a  path  constraint,  but  it  is  not  always  required, 
and  is  typically  problem  dependent. 

The  files  used  in  addition  to  the  main  program  file,  then  are: 

•  Cost  Function  -  outputs  are  the  endpoint  cost  and  running  cost.  In  the 
problems  being  examined,  only  endpoint  costs  are  required  and  are 
problem  specific.  In  the  following  cases  it  is  either  the  final  time, 
tf  (minimum  time  case),  or  the  final  mass,  mf  (minimum  fuel  case). 

•  Dynamics  Function  -  provides  the  differential  equation 
x(t)  =  f(x(t),u(t),t)  as  described  by  the  problem  dynamics.  The  same 
set  of  dynamics  equations  is  used  for  all  cases  studied. 

•  Events  Function  -  used  to  describe  the  boundary  conditions  to  include  the 
initial  and  final  states  of  the  problem.  For  example,  in  the  minimum- fuel, 
constrained-time  problem,  the  initial  state  is 
x0  =  [x°,  v0,z0,vv0,vv0,vz0,/«0,]  and  the  final  state  is 

xf  =  \jcJ  ,yf  ,zf  ,v/ ,vj  ,v/]  ,  noting  that  there  is  no  final  mass  specified 
since  that  is  the  parameter  being  solved  for. 

The  main  program  file  calls  these  files  and  then  requires  the  establishment  of 

upper  and  lower  bounds  on  the  states  and  the  controls.  This  is  important  step  for 

primarily  computational  reasons.  It  provides  a  range  of  values  that  the  DIDO©  algorithm 

can  use  in  the  output  of  state  and  control  trajectories.  Given  that  the  physical 

displacement  and  velocities  of  the  spacecraft  should  be  limited  to  a  region  near  the  earth 

and  moon,  the  upper  bounds  of  the  state  variables  should  be  provided  accordingly. 

However,  if  the  bounds  are  too  small,  the  resulting  trajectories  might  result  in  either  an 

infeasible  or  non-optimal  solution.  Excessively  large  bounds  will  start  to  affect 

computation  time.  DIDO©  also  allows  for  an  initial  guess  of  the  trajectories  that  gives 

the  program  a  starting  point,  however  it  is  possible  to  get  acceptable  results  without  an 
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initial  guess.  The  actual  algorithm  that  solves  the  problem  is  called  in  Matlab  by  a  single 
command  line  and  its  outputs  are  the  final  cost,  as  well  as  primal  and  dual  structures. 
The  primal  structure  includes  the  state  and  control  vectors  and  associated  discrete  times 
with  their  lengths  being  equal  to  the  number  of  LGL  nodes  used.  The  dual  structure 
includes  the  co-state  and  Hamiltonian  vectors,  also  equal  to  the  number  of  nodes.  These 
outputs  can  then  be  used  for  further  analysis  for  verification  of  optimality  and  feasibility. 

C.  BEST  PRACTICES 


1.  Scaling 


Most  computational  numerical  methods  make  use  of  unit  scaling,  without  loss  of 
generality,  to  increase  computational  efficiency  by  increasing  rate  of  convergence  while 
increasing  accuracy  (Ross  I.M.,  2009).  There  are  no  set  ground  rules  on  how  to  scale 
units,  however  a  unifonn  range  of  values  can  make  the  program  run  smoother  (Betts, 
2001).  In  astrodynamics  problems  in  particular,  scaling  becomes  an  important  factor  in 
problem  formulation  because  of  the  scale  at  which  spacecraft  are  operating.  Typically, 
the  range  of  values  for  the  distances  travelled  can  be  on  the  order  of  107  meters,  but  the 
timeframe  of  interest  may  only  be  on  the  order  of  a  few  hours.  Also  to  be  considered  are 
the  values  of  forces  in  kilonewtons,  for  example,  which  depend  on  the  standard  units  of 
mass,  distance,  and  time. 

For  the  problems  discussed  in  this  thesis,  the  following  “canonical”  units  were 
used  to  scale  the  variables  for  computation. 

•  1  Displacement  Unit  (DU)  =  1,737  km  (mean  radius  of  the  Moon) 


•  1  Mass  Unit  (MU)  =  20,339.9  kg  (wet  mass  of  spacecraft) 


DU 


g  l,033.9sec 


•  1  Time  Unit  (TU)  = -y  Mm  4,902.8W/sec" 

Where  juM  is  the  moon's  gravitational  constant. 


,  ..  ,  TT  1  DU  1, 737 km 

•  1  Velocity  Unit  = - = - =  1 .68  km  /  sec 

ITU  1,033.9  sec 


37 


1  Force  Unit  = 


\MU-\DU 

ITU2 


20,339.9tS.i,737^g33UA, 

(l,003.9sec) 


The  units  scaled  in  this  fashion  are  sufficiently  close  in  order  of  magnitude  to 
maintain  computational  efficiency  and  are  able  to  be  converted  to  and  from  the  actual 
(un-scaled)  values  for  analysis. 


2.  LGL  Nodes 


Typically,  these  types  of  complex  optimal  control  problems  cannot  be  solved  in  a 
single  step  with  a  large  number  of  LGL  nodes  primarily  due  to  computation  time.  The 
software  used  allows  the  user  to  choose  the  number  of  discrete  points  as  the  output  and 
allows  for  taking  the  output  of  a  lower  node  solution  and  using  it  as  a  guess  for  a  higher 
node  solution.  If  the  problems  are  well  scaled  and  the  problem  is  properly  formatted,  this 
allows  for  a  rapidly  converging  solution  without  excess  computational  time.  In  the 
problems  studied  in  this  thesis,  the  typical  starting  point  was  at  30  or  40  LGL  points  and 
then  each  incremental  step  increase  in  nodes  was  of  the  same  order.  The  exception  was 
the  study  of  the  high  thrust  singular  arc  in  which  a  much  lower  number  of  nodes  was  used 
to  start  because  of  the  short  time  duration  of  the  event. 


3.  Providing  an  Initial  Guess 

It  is  possible,  though  not  required,  to  provide  DIDO©  with  a  “guess”  trajectory 
for  the  controls  and  states.  Providing  a  reasonable  guess  will,  in  most  circumstances, 
produce  a  faster  result.  Each  of  the  cases  studied  in  this  thesis  started  with  a  low  number 
of  LGL  nodes  and  with  no  initial  guess,  however  as  the  number  of  nodes  and  subsequent 
fidelity  of  the  solutions  were  increased,  the  preceding  solutions  were  used  as  an  initial 
guess  for  the  next  iteration.  This  concept  of  using  an  initial  guess  was  used  in  two 
different  ways  depending  on  the  situation. 

Starting  at  the  low  fidelity  solution,  typically  30  to  40  LGL  nodes,  it  is  possible  to 
take  the  resulting  primal  structure  (controls,  discrete  time  values,  and  states)  and  use 
them  as  the  initial  guess  for  a  subsequent  run  of  DIDO©  at  a  higher  number  of  nodes. 

Experience  shows  that  increasing  in  increments  of  no  more  that  30  to  40  works  best. 
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This  method  is  called  bootstrapping,  but  cannot  be  used  without  analysis  of  the  preceding 
result  so  as  not  to  propagate  errors  caused  by  bad  problem  formulation,  scaling,  etc. 
Therefore,  each  step  of  the  problem  solution  process  should  be  checked  for  optimality 
conditions,  feasibility,  and  other  verification  and  validation  techniques  discussed  below. 

The  difficulty  with  the  angles  is  that  as  previously  mentioned  they  are  unbounded 
controls,  however,  in  coding  the  problem  it  is  typically  unwise  to  allow  infinite  bounds 
for  computational  reasons.  The  problem  as  coded  must  establish  small  enough  bounds  so 
that  it  can  be  handled  by  the  computer,  but  must  also  be  large  enough  to  facilitate  the 
possible  range  of  solutions.  Witha,/?  e  [— 2tt,  2tt]  ,  this  allowed  for  a  wide  range  of 
angles  from  which  the  solution  could  be  found. 

The  resulting  control  angles  are  allowed  within  the  program  to  reach  the  limits  of 
the  bounded  region,  however  when  this  happens  it  is  difficult  to  detennining  whether  the 
control  has  hit  a  “physical”  limit  imposed  by  the  code  or  the  angle  is  equal  to  that  limit, 
which  in  this  case  is  zero  degrees.  Because  of  the  non-unique  nature  of  trigonometric 
functions,  a  one-to-one  mapping  of  solutions  is  not  possible  and  computer  algorithms 
typically  cannot  discern  the  difference  between  them  (i.e.  sin(;r/2)  =  sin(5;r/2)  =  1 ). 
Therefore,  for  any  computer  algorithm  additional  conditions  may  have  to  be  imposed. 

In  certain  cases,  angle  reduction  formulas  could  be  used  to  convert  the  output 
angles  into  a  set  of  angles  that  fall  between  a,fie  \-n,n\  which  would  then  provide  the 

full  range  necessary  for  the  problem  solution.  From  Equation  (1.16),  the  Cartesian 
components  of  the  thrust  vector  can  be  re-written  as: 


Tx=T  cos(a)cos(/?) 

Ty  =  Esin  («)  cos  (/?)  (4.1) 

Tz=Tsm(p) 


From  these  equations,  explicit  formulas  for  the  angles  can  be  derived. 
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a  =  tan 


( 

V 


T  'j 

y 

TJ 


ft  =  sin  1 


71  = 


TJ  +T2  +t2 


(4.2) 


Noting  that  a  singular  condition  occurs  when  F  =  0 ,  this  can  be  neglected  by  the 

fact  that  any  combination  of  angles  has  no  effect  on  a  spacecraft’s  trajectory  when  the 
thrust  is  zero,  therefore  when  this  is  the  case  the  angles  can  be  forced  to  zero  with  no  loss 
of  generality  to  the  problem.  Figures  14  and  15  show  how  a  set  of  resulting  angles  from 
DIDO©  can  be  recalculated  using  Equation  (4.2)  which  can  then  be  used  as  the  initial 
guess  for  the  angles  for  the  next  iteration. 


Figure  14.  Recalculation  of  Azimuth  Control  Angle 
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8 


Time  (hrs) 

Figure  15.  Recalculation  of  Elevation  Control  Angle 

As  the  problem  increases  in  fidelity  and  the  control  angles  start  to  fall  within  a 
tighter  range  of  values,  it  is  then  possible  to  tighten  up  the  coded  bounds  to  help  speed  up 
computation. 

D.  VERIFICATION  AND  VALIDATION 
1.  Examining  Initial  Output 

One  of  the  first  steps  in  the  analysis  conducted  was  to  generate  plots  of  the  state 
and  control  trajectories  as  well  as  the  evolution  of  the  Hamiltonian  and  the  co-states.  One 
of  the  first  checks  is  to  determine  if  any  unexpected  or  undesired  limits  of  the  imposed 
bounds  were  met.  In  the  case  of  thrust,  the  limits  are  expected  because  of  the  actual 
physical  limitations  of  the  thruster  (no  negative  thrust  and  a  maximum  attainable  thrust). 
Hitting  the  bounds  on  the  control  angles  was  covered  above.  For  a  minimum  time 
problem,  the  final  time  should  not  hit  the  upper  bound  because  it  usually  results  in  an 
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infeasible  solution  indicating  that  the  time  horizon  should  be  longer.  In  the  minimum 
fuel  problem,  the  upper  time  bound  is  typically  fixed,  and  therefore  coded  appropriately. 

2.  Checking  Necessary  Conditions  for  Optimality 

From  the  dual  variables  output  of  DIDO©  the  Hamiltonian  and  co-state  vectors 
can  be  plotted  with  respect  to  time  and  using  the  necessary  conditions  from  Chapter  III 
can  be  used  to  verify  the  optimality  of  the  solution. 

3.  Checking  Feasibility 

The  DIDO©  endpoint  conditions  can  be  validated  against  the  given  endpoint 
conditions,  however  the  purpose  of  perfonning  the  feasibility  analysis  is  to  validate  the 
DIDO©  solution  and  detennine  if  the  solution  has  enough  fidelity  to  be  “flyable.”  A 
method  of  verifying  the  feasibility  of  the  solution  is  to  take  the  same  initial  conditions  for 
the  problem  and  propagate  the  same  dynamics,  using  Matlab’s  ode45  solver,  over  the 
same  timeframe  to  evaluate  how  closely  the  state  trajectories  follow  the  DIDO©  solution. 
This  requires  an  interpolation  of  the  generated  controls  for  areas  between  the  nodes  since 
there  many  more  points  used  in  the  ode45  solver  compared  to  the  sparse  number  of  points 
generated  by  DIDO©.  Plotting  the  DIDO©  and  propagated  state  trajectories 
simultaneously  can  provide  a  good  visual  indicator  of  whether  the  solution  is  making 
physical  sense,  assuming  the  dynamics  are  correct.  However,  for  problems  that  cover  a 
large  time  horizon,  the  results  will  tend  to  diverge  as  time  goes  on  and  this  is  primarily  a 
result  of  interpolation  and  propagation  errors  inherent  to  numerical  solutions. 

A  technique  that  has  been  proven  to  be  successful  in  reducing  the  propagation 
errors  of  the  interpolated  controls  from  the  optimal  solution  generated  by  DIDO©  is 
called  the  Bellman  Pseudospectral  (PS)  Method  (Ross,  I.M.,  Gong,  Q.,  and  Sekhavat,  P., 
2008).  Its  application  is  based  on  the  Bellman  Principle  of  Optimality  that  states  that  for 
an  optimal  path  between  two  points,  any  point  taken  along  the  path  will  follow  the 
optimal  path  to  the  final  point. 

The  Bellman  Principle  of  Optimality  states  that  along  a  given  optimal  trajectory,  a 

segment  starting  from  some  intennediate  point  on  that  trajectory  and  ending  at  the 
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original  termination  point  will  also  be  optimal.  If  some  trajectory  AB  meets  the 
necessary  conditions  for  optimality,  then  a  segment  starting  from  point  C  on  the  path  will 
follow  the  same  optimal  path  to  the  final  point.  In  Figure  16  the  total  cost  of  the  optimal 
trajectory  AB  is  J0 .  If  the  point  dies  on  the  trajectory  A3  ,  then  the  trajectory  from  CB 

will  also  be  an  optimal  path  on  AB  ,  and  J0  =  Jx  +  J2  (Kirk,  1970). 


A 


Figure  16.  Bellman  Principle 

Yan  et  al.  use  the  Bellman  PS  Method  as  a  means  for  mesh  refinement  and  taking 
advantage  of  DIDO©  node  distribution  which  puts  a  higher  concentration  of  nodes  near 
the  beginning  and  end  of  the  time  horizon  of  the  optimal  trajectories.  Hence,  the 
resolution  of  data  near  the  beginning  and  end  of  the  trajectory  is  higher  than  in  areas  in 
the  middle.  Increasing  the  number  of  nodes  would  only  offer  a  limited  increase  in 
resolution  near  the  middle  of  the  trajectory  which  would  come  at  a  high  cost  in 
computational  efficiency.  They  found  that  by  applying  the  technique  in  this  manner,  the 
interpolation  and  propagation  errors  were  significantly  mitigated  by  a  few  orders  of 
magnitude  (Yan  et  al.,  2010).  The  application  of  the  corollary  statement  was  used  in  the 
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detailed  analysis  of  the  singular  arc  which  used  two  points  on  the  trajectory  prior  to  and 
after  the  second  TEI  maneuver  of  the  main  engine  solution  as  a  segment  which  lay  on  an 
existing  optimal  path. 

4.  Comparison  to  Previous  Studies 

Finally,  the  results  of  the  case  studies  in  this  thesis  can  be  compared  to  the  results 
by  Yan  et  al.  study  of  optimal  moon  to  earth  trajectories  using  the  Orion  spacecraft  main 
engines  as  well  as  the  study  conducted  by  Park  et  al.  on  numerical  verification  of  singular 
arcs  for  the  same  problem  using  a  low  thrust  engine.  The  rocket  engine  parameters, 
resulting  fuel  consumption,  and  control  maneuvers  can  be  compared  for  analysis. 
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V.  AUXILIARY  ENGINE  FOR  MOON  TO  EARTH  RETURN 

MISSION 


For  the  three  cases  for  the  moon  to  earth  return  mission,  the  scenario  depicted  is  a 
loss  of  the  main  engine  with  the  same  mission  time  requirement  and  fuel  constraint.  The 
initial  position  and  velocity,  or  state,  of  the  spacecraft  is  at  a  designated  point  in  a  low- 
lunar  orbit  from  which  the  mission  commences.  The  final  state  achieves  the  conditions 
that  will  take  the  spacecraft  to  the  earth  interface  point  without  any  further  burn 
maneuvers.  The  initial  epoch  is  given  as  April  4,  2024,  15:30:00  TDT.  Each  of  the  three 
cases  studied,  the  final  position  and  velocity  remain  the  same  and  the  results  of  changing 
the  boundary  conditions  on  time  and  fuel  consumed  are  examined  for  the  purpose  of 
determining  the  overall  feasibility  of  using  the  auxiliary  engines  as  an  alternative  to  the 
main  engines  for  the  return  mission. 

The  spacecraft  characteristics  and  fixed  initial  and  tenninal  conditions  are  as 
follows  (Scarritt,  S.,  Marchand,  B.,  and  Weeks,  M.,  2009): 

•  Initial  Mass:  20,339.9  kg  (wet  mass) 

•  Total  Fuel  Mass:  8,063.65  kg 

•  Auxiliary  Engine  Thrust:  4,448.0  N 

•  Auxiliary  Engine  Isp:  309  sec 

•  Initial  State  in  the  J2000  Moon  centered  inertial  frame: 


x0 

-1,236.7970783385588  km 

To 

1,268.1142350088496  km 

468.38317094160635  km 

Vx0 

0.0329108058365355  km/ sec 

v,o 

0.589269803607714  km  /sec 

_V--oJ 

-1.528058717568413  km! sec 

45 


Final  State  (to  reach  specified  earth  interface  point) 


xf 

2,106.40349727063  km 

yf 

3,596.40440859595  km 

zf 

608.167937792834  km 

3 

-0.151018343429892  bn /sec 

Vyf 

-0.321178999121464  km /sec 

3/. 

_  -1.79844521633271  km /sec_ 

A.  THE  MINIMUM  TIME,  FUEL  FREE  PROBLEM 

1.  Problem  Formulation 

The  first  case  of  this  problem  was  formulated  as  a  minimum  time  solution  with 
the  final  mass  allowed  to  be  free,  but  limited  to  a  lower  bound  of  10%  of  the  total 
spacecraft  wet  mass.  While  this  problem  did  not  satisfy  the  fuel  requirements,  it  was  an 
important  first  step  in  determining  an  overall  feasible  engineering  solution.  By  making 
both  the  time  and  fuel  constraints  ‘free,’  it  was  possible  to  first  detennine  the  feasibility 
of  the  mission  time  requirement  imposed  on  this  segment  of  the  return  mission.  With 
time  as  the  cost  function  to  be  minimized,  the  spacecraft  was  allowed  to  bum  as  much 
fuel  as  needed,  even  beyond  the  maximum  fuel  capacity  of  8,063.65  kg.  In  the  Matlab 
code,  however,  an  artificial  lower  limit  of  10%  of  the  spacecraft  dry  mass  was  imposed  to 
prevent  the  mass  from  going  negative.  Provided  an  optimal  solution  existed,  the  goal  was 
to  determine  if  the  minimal  final  time  was  less  than  the  48-hour  mission  constraint.  This 
in  turn  implied  that  the  use  of  the  specified  auxiliary  engine  would  be  feasible,  at  least  in 
tenns  of  the  mission  time  constraint. 
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Minimum  Time  Formulation 


Where  x  = 


x,  y,  z, 


Vx’  V 


i  and  u=\T, a, /?], 


Minimize 
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rs  m 

vsiiziA  +  T 

rs3  m 
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rs3  m 


2.  Results  and  Analysis 

This  case  was  solved,  starting  at  30  LGL  nodes  and  bootstrapping  the  solution  up 
to  80  LGL  nodes  in  order  to  achieve  a  smooth  and  convergent  solution.  As  expected,  in 
order  to  achieve  the  desired  end  point  condition  in  minimum  time,  the  spacecraft  will 
bum  fuel  at  a  maximum  rate,  using  maximum  thrust,  throughout  the  maneuver  as 
depicted  in  Figure  17.  Figures  18  and  19  show  that  the  azimuth  and  elevation  control 
angles  are  continuous  and  while  they  are  not  physical  constraints  on  the  problem,  the 
angles  do  not  hit  the  bounds  specified  in  the  Matlab  code.  While  the  fuel  mass  consumed 
exceeds  the  actual  fuel  capacity  of  the  spacecraft,  the  total  time  of  the  maneuver  is  much 
less  than  the  48  hours  required  by  the  mission.  This  implies  that  the  auxiliary  engines, 
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subject  to  the  dynamics  of  the  problem,  can  be  used  and  will  satisfy  the  time 
requirements  of  the  mission.  Further  analysis  follows  to  verify  that  the  auxiliary  engines 
can  be  used  to  meet  the  fuel  requirements  as  well. 
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Figure  17.  Thrust  Profde 
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Figure  18.  Azimuth  Control  Profile 


Figure  19.  Elevation  Control  Profile 
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a.  Feasibility  of  the  Solution 

As  part  of  the  feasibility  analysis,  the  solution  is  propagated  using  the 
ode45  function  in  Matlab,  using  the  initial  boundary  condition  and  interpolating  the 
controls  over  the  course  of  the  trajectory  and  comparing  the  continuous  solution  against 
the  DIDO©  solution  of  discrete  points.  By  examining  the  profiles  of  the  state  vectors  for 
position,  velocity,  and  mass,  it  was  found  that  the  DIDO©  solution  follows  the 
propagated  solution  and  is  deemed  a  feasible  solution.  Figures  20,  21,  and  22  show  the 
displacement,  velocity,  and  mass  trajectories.  Table  3  shows  the  differences  between  the 
DIDO©  and  the  propagated  solutions.  While  the  DIDO©  solution  achieves  the  specified 
boundary  conditions,  the  propagated  solution  shows  an  error  greater  than  the  required 
limits  of  1  km  for  displacement  and  1  m/s  for  velocity. 


Figure  20.  Displacement  Vector 


50 


2.5 


A  _ I _ I _ I _ I _ I _ 

0  0.5  1  1.5  2  2.5  3 

Time  (hrs) 


F igure  21.  T angential  Velocity  V ector 


Figure  22.  Fuel  Consumption 
51 


Final  Position  Error  (km) 

Final  Velocity  Error  (m/s) 

DIDO© 

4.5e-13 

2.2e-13 

Propagated  Solution 

17.0 

4.2 

Table  3.  DIDO  vs.  Propagated  Solution 


The  discrepancies  in  the  final  position  and  velocities  are  due  to  the 
numerical  sensitivities  of  this  type  of  problem.  Very  slight  interpolation  errors  at  the  start 
of  the  orbit  transfer  can  cause  very  large  errors  at  the  terminal  position  (Yan  et  ah,  2010). 
In  this  particular  problem,  the  total  time  horizon  is  shorter  than  the  48-hour  requirement 
and  thus  results  in  a  relatively  small  error.  Within  two  iterations  of  the  Bellman 
technique,  the  required  end-point  errors  are  achieved.  For  brevity,  the  Bellman  technique 
is  demonstrated  for  this  problem  only,  but  can  be  used  to  resolve  the  difference  between 
the  DIDO©  solution  and  the  propagated  solution  in  the  minimum  fuel  case. 

The  Bellman  pseudo-spectral  (PS)  technique  is  summarized  as  follows 
(Ross  et  ah,  2008): 

1 .  Solve  the  optimal  control  problem  using  a  reasonably  low  number 
of  LGL  nodes  to  generate  a  discrete  time  solution. 

2.  Partition  the  time  interval  j^ty^into  Bellman  segments 

t°  <  tl  < ...  <  t v"  =  t" .  The  segments  do  not  need  to  be  uniformly 
spaced,  however  extend  from  a  specified  time  until  the  terminal 
time,  tf . 

3.  Propagate  the  system  dynamics  from  t° to  fusing  T0 as  the  initial 
condition  and  any  method  of  continuous-time  reconstruction  of  the 
controls,  wj  (t),t  e  [j0,?1]  . 

4.  Set  x0  =  x1  (V )  and  t°  =  tl  and  go  back  to  Step  1 . 

This  algorithm  ends  when  the  final  required  conditions  are  met. 
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5. 


Using  the  Bellman  technique  in  two  iterations,  the  final  required  accuracy 
was  attained.  The  first  Bellman  segment  ^,,^Jwas  set  to  start  at  0.8363  hours.  For 
simplicity,  a  time  associated  with  one  of  the  discrete  time  points  was  chosen.  For  the 
second  Bellman  segment  \  t2,tf  J ,  a  start  time  of  0.9640  hours  was  chosen.  Table  4  shows 

the  segment  partition  and  the  resulting  final  position  and  velocity  errors.  After  two 
iterations,  the  required  errors  for  both  position  and  velocity  are  attained. 


Position  Error  (km) 

Velocity  Error  (m/s) 

Without  Bellman  technique 

17.0 

4.2 

First  Bellman  segment 

1.2 

0.31 

Second  Bellman  segment 

0.29 

0.16 

Table  4.  Effects  of  the  Bellman  PS  Technique 


Another  verification  of  the  feasibility  of  the  solution  is  to  examine  the 
osculating  orbital  elements  throughout  the  maneuver.  Figure  23  shows  the  evolution  of 
the  changing  semi-major  axis  throughout  the  maneuver.  Towards  the  end  of  the 
maneuver,  the  semi-major  axis  goes  to  infinity  which  is  consistent  with  a  parabolic 
trajectory,  but  then  completes  the  maneuver  at  a  finite  semi-major  axis  value.  Figure  24 
shows  the  eccentricity  constantly  changing  throughout  the  orbit  and  eventually  passes 
through  e  =  1,  a  parabola,  and  ends  on  a  value  greater  than  one  which  defines  a 
hyperbolic  trajectory.  Both  parabolic  and  hyperbolic  trajectories  imply  that  a  spacecraft 
has  sufficient  energy  to  overcome  the  gravity  well  of  the  body  which  it  has  been  orbiting 
(Curtis,  2005).  This  result  shows  that  the  spacecraft  in  this  case  has  escaped  the  moon’s 
gravity.  Finally,  Figure  25  shows  a  continual  change  in  the  inclination  of  the  trajectory 
with  respect  to  the  moon’s  equator.  The  final  states  of  the  semi-major  axis,  eccentricity, 
and  inclination  are  consistent  with  the  end  point  conditions  found  for  the  minimum  fuel 
solution  for  the  main  engines  (Yan  et  ah,  2010). 
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Figure  23.  Semi-major  Axis 


Figure  24.  Eccentricity 
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Figure  25.  Inclination 


b.  Verification  of  the  Optimality  of  the  Solution 

The  necessary  conditions  for  optimality  are  met  and  are  examined  as 
prescribed  in  Chapter  III.  In  Figure  26,  the  Hamiltonian  is  constant  throughout  the 
maneuver  and  its  value  is  consistent  for  a  minimum  time  solution.  From  the 
transversality  condition,  the  co-state  related  to  mass  has  a  terminal  value  equal  to  zero  as 
shown  in  Figure  27. 
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Figure  26.  Flamiltonian  Evolution 


Figure  27.  Co-state  (Xm)  Evolution 
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The  condition  derived  from  the  HMC  that  the  angle  of  thrust  is  opposite 
the  velocity  co-state  is  depicted  in  Figure  28.  During  the  maneuver,  the  angle  maintains  a 
value  of  180  degrees  as  the  thrust  is  “on”  throughout.  And  finally,  the  Switching 
Function  as  defined  by  the  HMC  should  be  negative  if  the  constrained  control,  thrust,  is 
at  maximum  which  is  true  in  this  case.  Figure  29  shows  that  the  switching  function  is 
indeed  negative  throughout  the  maneuver  which  is  consistent  with  the  thrust  being  at  its 
maximum  value  throughout  the  maneuver. 


Figure  28.  Primer  Angle 
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Figure  29.  Switching  Function 


3.  Conclusion  for  Minimum  Time  Case 

The  results  show  that  the  solution  determined  in  the  minimum  time  problem  is 
optimal  with  respect  to  time.  The  maneuver  achieves  the  final  boundary  conditions  well 
within  the  48-hour  time  constraint  for  the  TEI  maneuvers  for  the  moon  to  earth  return 
mission,  which  informs  that  the  auxiliary  engines  are  capable  for  meeting  the  mission 
time  requirements.  The  next  step  is  to  determine  the  minimum  fuel  required  that  satisfies 
the  mission  time  requirements  and  determine  the  feasibility  of  using  the  auxiliary  engines 
subject  to  the  fuel  capacity  constraints  of  the  spacecraft. 
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B.  THE  MINIMUM  FUEL,  FIXED  TERMINAL  TIME  PROBLEM 

1.  Problem  Formulation 

The  auxiliary  engine  specifications  as  well  as  the  initial  and  terminal  boundary 
conditions  remain  the  same  as  in  the  minimum  time  optimal  control  problem  formulation. 
In  this  case,  the  problem  is  formulated  as  to  minimize  the  fuel  consumption  over  a  fixed 
time  horizon.  The  time  horizon  is  constrained  to  48  hours  and  represents  the  allowable 
window  for  the  Orion  spacecraft  to  conduct  the  TEI  burns  for  an  expedited  return  to 
earth.  The  goal  was  to  find  a  feasible  solution  to  the  problem  of  conducting  the  return 
mission  using  solely  the  auxiliary  engines.  Because  of  the  rocket  characteristics,  the 
auxiliary  engine  uses  more  fuel  for  a  given  duration  of  burn  than  does  the  main  engine. 
This  solution,  if  found,  was  to  show  that  the  spacecraft  would  have  enough  fuel  onboard 
to  conduct  the  return  mission.  Additionally,  verification  of  the  optimality  of  the  solution 
is  necessary. 


59 


Minimum  Time  Formulation 
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2.  Results  and  Analysis 

Figures  30,  31,  and  32  show  the  resulting  control  profiles  for  the  minimum  fuel 
problem.  The  problem  was  initially  started  using  30  LGL  nodes  and  then  bootstrapping 
the  solution  to  a  final  solution  using  300  LGL  nodes.  The  goal  was  to  balance 
computation  time  by  limiting  the  number  of  nodes  at  the  same  time  as  achieving  a 
smooth  and  viable  control  solution.  As  the  number  of  LGL  nodes  was  increased  for  each 
iteration  of  the  solution,  the  control  angles  were  forced  to  zero  whenever  the  thrust  was 
zero.  This  helped  speed  up  the  computation  of  the  DIDO©  results  and  was  assumed  that 
this  would  have  no  effect  on  the  optimal  solution  since  the  control  angles  only  influence 
the  spacecraft  trajectory  whenever  there  is  a  thrust  associated  with  them.  The  thrust 
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Thrust  (KN 


profile  is  similar  to  the  main  engine  trajectory  solved  by  Yan  et  al.  in  that  it  contains  three 
distinct  TEI  maneuvers,  with  the  second  burn  having  the  characteristics  of  a  singular  arc. 


Figure  30.  Thrust  Control  Structure 
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Figure  3 1 .  Azimuth  Control  Angle 


Figure  32.  Elevation  Control  Angle 
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The  entire  series  of  maneuvers  are  completed  within  the  bounded  time  horizon 
and,  as  depicted  in  Figure  33,  uses  a  total  of  7,215.6  kg  of  fuel,  less  than  the  given  fuel 
capacity  of  the  spacecraft  which  is  8,063.65  kg. 


Figure  33.  Fuel  Consumption 


a.  Feasibility  of  the  Solution 

A  feasibility  analysis  of  the  solution  was  conducted,  similarly  to  the 
previous  case,  where  the  solution  is  propagated  using  the  Matlab  ode45  function.  The 
controls  are  interpolated  over  the  propagated  dynamics  to  see  how  well  the  discrete 
DIDO©  solution  maps  to  the  integrated  solution.  Figures  34  and  35  illustrate  that  the 
displacement  and  tangential  velocities  of  the  discrete  solution  follow  the  propagated  state 
trajectories.  Some  divergence  exists  towards  the  end  of  the  trajectory  due  to  numerical 
propagation  errors.  The  difference  in  the  terminal  condition  of  the  propagated  solution 
and  the  DIDO©  solutions  are  502. 1  km  in  displacement  and  0. 10  km/sec  in  velocity.  The 
variation  in  final  mass  was  approximately  0.1  kg.  While  DIDO©  finds  a  solution  that 
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meets  the  given  boundary  conditions,  a  resolution  to  the  validation  and  verification  of  the 
results  need  be  conducted  since  the  propagated  solution  does  not  meet  the  displacement 
error  and  velocity  error  requirements  of  1  km  and  1  m/s  respectively.  These  errors  can  be 
eliminated  using  the  Bellman  technique  described  in  the  minimum  time  problem. 
However,  the  initial  analysis  shows  that  since  the  propagated  solution  and  discrete 
solutions  are  closely  matched,  the  solution  appears  feasible. 


Figure  34.  Displacement  State  Trajectory 
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Figure  35.  Velocity  State  Trajectory 

Again,  like  the  analysis  perfonned  in  the  minimum  time  problem,  the 
evolution  of  the  osculating  orbital  elements  was  investigated.  In  Figure  36,  the  semi¬ 
major  axis  changes  primarily  due  to  the  first  and  last  TEI  burns.  The  first  burn  is  an 
apoapsis  raising  maneuver  and  the  final  burn  provides  the  kinetic  energy  required  for  the 
spacecraft  to  achieve  escape  velocity  from  the  moon.  In  a  parabolic  trajectory,  the  semi¬ 
major  axis  is  not  defined,  hence  the  singular  spike  in  the  figure.  Figure  37  shows  that  the 
eccentricity  of  the  orbit  increases  to  a  value  greater  than  one,  which  means  that  the 
spacecraft  departs  the  moon’s  orbit  on  a  hyperbolic  trajectory,  requiring  a  transition 
through  a  parabolic  trajectory  where  the  eccentricity  equals  exactly  one.  Figure  38  shows 
an  inclination  change  of  about  15  degrees  corresponding  to  the  second,  singular  arc 
maneuver.  The  terminal  values  of  the  osculating  orbital  elements  match  those  of  the 
main  engine  TEI  maneuver  sequence  (Yan  et  ah,  2010).  It  should  be  noted  that  in  each  of 
the  TEI  bums,  a  combination  of  inclination  and  eccentricity  changes  is  found.  This 
appears  to  be  the  case  since  each  of  the  maneuvers  is  finite  in  duration,  rather  than 
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impulsive.  This  effect  is  amplified  because  of  the  lower  thrust  engine  which  requires  a 
longer  burn  time  to  perform  a  given  maneuver  than  one  with  higher  thrust. 
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Figure  36.  Semi-major  Axis  With  Respect  to  Moon  Centered  Frame 
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Figure  37.  Eccentricity  With  Respect  to  Moon  Centered  Frame 
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Figure  38.  Inclination  With  Respect  to  Moon  Centered  Frame 
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A  final  check  to  ascertain  the  feasibility  of  the  solution  is  to  compare  the 
final  fuel  consumption  to  the  results  from  several  authors  who  solved  the  same  minimum 
fuel  earth  return  problem  using  the  same  initial  and  final  boundary  conditions,  but  used 
different  engine  characteristics.  In  Yan  et  ah,  as  previously  described  solve  the  problem 
using  the  main  engine  with  a  thrust  of  33.6  kN  and  an  Isp  of  326  seconds.  Park  et  al. 
used  a  rocket  engine  with  a  thrust  of  3  kN  and  an  Isp  of  326  seconds.  The  results  in 
Table  5  show  the  resulting  fuel  consumption  of  the  present  problem  using  only  the 
auxiliary  engines  for  the  entire  return  mission  compared  to  each  of  these  cases.  The  fuel 
consumption  value  of  7,215.6  kg  should  be  close  to  the  value  reached  by  Park  with 
similar  engine  characteristics  and  greater  than  the  main  engine  which  with  a  higher  thrust 
can  perfonn  the  maneuvers  with  shorter  burn  durations  and  greater  fuel  efficiency 
(greater  Isp). 


Engine  Configuration 

Isp  (sec) 

LGL  Nodes 

Fuel  Consumed  (kg) 

Auxiliary  Engines  (4.4  kN) 

309 

300 

7,215.6 

Auxiliary  Engines  (3  kN) 

326 

160 

7,245.1 

Main  Engines  (33.6  kN) 

326 

500 

6,657.3 

Table  5.  Fuel  Consumption  Comparison 


b.  Verification  of  the  Optimality  of  the  Solution 

Figures  39  and  40  show  the  Hamiltonian  constant  at  a  value  of  zero,  and 
the  transversality  condition  of  the  co-state  related  to  mass  equal  to  -1  at  the  terminal 
point. 
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Figure  39.  Flamiltonian  Evolution 


Figure  40.  Co-state  (Xm)  Evolution 
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Figure  41  shows  the  evolution  of  the  switching  function  over  the  entire 
maneuver  sequence.  The  switching  function  is  negative  for  maximum  thrust  and  positive 
for  zero  thrust,  however  at  the  singular  arc  the  switching  function  also  appears  to  be 
negative.  The  singular  arc  is  demonstrated  to  be  optimal  by  examining  the  derivatives  of 
the  switching  function  up  to  the  fourth  derivative.  Figure  42  shows  that  the  fourth 
derivative  is  equal  to  zero  during  the  singular  arc.  The  generalized  Legendre-Clebesch 
convexity  condition  during  the  singular  arc  burn  is  shown  in  Figure  43.  Figure  44  is 
included  solely  to  demonstrate  that  it  is  possible  to  have  regions  where  the  Legendre- 
Clebesch  condition  is  non-positive. 
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Figure  4 1 .  Thrust  and  Switching  Function 
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Generalized  Legendre-Clebesch  (scaled) 


Time  (hrs)  Time  (hrs) 


Figure  42.  Derivatives  of  the  Switching  Function  During  Singular  Arc 


Figure  43.  Generalized  Legendre-Clebesch  Condition  During  Singular  Arc 
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Figure  44.  Generalized  Legendre-Clebesch  Condition  Over  Entire  Trajectory 


xIO"3 


?"  ^^Negative  Value 
• 

Figures  45,  46,  and  47  show  in  detail  the  features  of  the  TEI  bums;  they 
show  the  thrust  and  control  angles  and  also  that  for  each  burn,  the  primer  angle  between 
the  thrust  and  the  velocity  co-vector  is  equal  to  180  degrees  as  necessary  to  meet  the 
Hamiltonian  Minimization  Condition.  With  the  exception  of  the  sparse  nature  of  the 
control  trajectory  of  the  second,  singular  TEI  bum,  the  control  solution  appears  smooth 
enough  for  a  control  solution.  From  a  practical  control  perspective,  however,  it  is  desired 
to  increase  the  fidelity  of  the  discrete  solutions  into  a  more  continuous  solution. 
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Figure  45.  First  TEI  Burn 


Figure  46.  Second  TEI  Burn 
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Figure  47.  Third  TEI  Bum 


As  pointed  out  in  the  feasibility  analysis,  the  final  displacement  and 
velocity  errors  are  due  to  the  interpolation  errors  propagated  over  a  long  time  horizon. 
Figures  48  and  49  show  how  the  error  between  the  DIDO©  solution  and  the  interpolated 
solution  increases  over  the  entire  solution  trajectory.  Using  the  Pseudospectral  Bellman 
Technique  used  by  Ross  et  ah,  and  demonstrated  in  the  minimum  time  problem,  it  is 
possible  to  minimize  the  errors  caused  by  interpolation  of  the  controls. 
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Figure  48.  Displacement  Error 


Figure  49.  Velocity  Error 
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c. 


CONCLUSIONS  FOR  USE  OF  AUXILIARY  ENGINES  FOR  RETURN 
MISSION 


The  resulting  optimal  control  trajectory  shows  that  in  the  event  of  a  total  main 
engine  failure,  the  auxiliary  engines  could  be  used  for  the  entire  return  mission,  meeting 
the  48-hour  minimum  time  requirement  to  complete  the  third  TEI  maneuver. 
Additionally,  the  fuel  optimal  solution  shows  that  7,215.6  kg  is  the  theoretical  minimum 
required  fuel  to  be  carried  onboard,  not  including  any  design  margins.  This  result  then 
potentially  reduces  the  total  amount  of  fuel  required,  thus  freeing  up  additional  mass  for 
other  aspects  of  the  spacecraft  design. 
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VI.  MAIN  ENGINE  SINGULAR  ARC  ANALYSIS 


The  Yan  et  al.  results  for  the  minimum  fuel  solution  for  the  earth  return  mission 
using  the  main  engines,  they  found  a  singular  arc  as  part  of  the  control  trajectory  for  the 
second  TEI  maneuver.  In  their  solution,  the  singular  bum  is  made  up  of  only  three 
distinct,  non-zero  thrust  points.  Since  a  high  fidelity  control  solution  is  desired,  it  would 
be  advantageous  to  have  additional  points  from  which  to  approximate  a  continuous 
solution.  As  a  secondary  investigation,  NASA  is  interested  in  detennining  the  feasibility 
of  using  the  auxiliary  engines  in  place  of  the  main  engines  for  the  singular  arc  TEI 
maneuver.  Therefore,  the  purpose  of  the  singular  arc  study  is  to: 


(1)  Generate  a  higher  fidelity  control  solution  from  the  existing  singular  arc 
data  generated  by  Yan  et  al. 

(2)  Determine  the  feasibility  of  using  the  auxiliary  engines  for  the  singular  arc 
and  ascertain  the  fuel  penalty  if  any 

(3)  Verily  and  validate  the  necessary  conditions  for  optimality  for  any  feasible 
solution(s)  generated 


The  separation  between  the  nodes  from  the  original  main  engine  solution  is  due 
primarily  to  its  location  along  the  entire  trajectory  solution.  The  output  data  points  from 
DIDO©  are  non-unifonnly  spaced,  being  more  dense  at  the  beginning  and  end  of  the 
trajectory  and  sparse  in  the  middle.  The  singular  arc  was  in  the  sparse  region  of  the 
original  solution  covering  a  time  span  of  approximately  29  minutes.  The  data  points 
around  the  singular  arc  were  spaced  approximately  7  minutes  apart  as  compared  to  the 
points  around  the  final  TEI  bum  had  a  spacing  range  of  approximately  0.4  to  0.8  minutes 
between  them.  The  approach  therefore  was  to  examine  the  singular  arc  by  fixing 
boundary  conditions  close  to  and  around  the  start  and  completion  of  the  maneuver.  Since 
the  goal  is  to  increase  the  number  of  discrete  points  for  a  solution  with  higher  resolution, 
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the  problem  uses  a  small  number  of  LGL  nodes  to  compute  a  solution  and  then  using  the 
bootstrap  technique  to  increase  the  nodes  until  the  control  trajectories  converge  in  a 
smooth  fashion. 

Recalling  that  the  Bellman  Principle  of  Optimality  states  that  along  a  given 
optimal  trajectory,  a  segment  starting  from  some  intermediate  point  on  that  trajectory  and 
ending  at  the  original  termination  point  will  also  be  optimal.  It  says  nothing  about  an 
intennediate  segment.  Suppose  there  are  two  points,  C  and  D ,  forming  a  segment  along 
the  same  optimal  trajectory.  Even  though  the  original  optimal  trajectory  can  pass  through 
these  points,  it  does  not  mean  that  the  original  path  on  AB  between  them  is  optimal  on 
the  path  between  C  and  D  .  There  exists  an  optimal  path  CD  that  may  not  lie  on  AB  as 
shown  in  Figure  50.  If  the  minimized  cost  along  AB  is  J(] ,  and  the  sum  of  the  segments 
on  AB  is  J0  =  Jl  +  J2  +  J:, ,  then  the  cost  of  an  optimal  trajectory  CD  not  on  AB  must  be 
J2  >  J2 .  Put  another  way,  Jx  +  J2  +  J,  <  ./,  +J'2+  . 


A 


Figure  50.  Bellman  Curve  and  Segment 
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The  singular  arc  generated  from  the  Yan  et  al.  solution  can  be  said  to  lie  on  a 
segment  of  a  fuel-optimal  trajectory.  In  order  to  maintain  the  same  cost  along  that 

segment,  so  that  from  the  example  above,  J2  =  J2 ,  then  the  dashed  line  trajectory  must  lie 

on  the  solid  line  CD .  This  implies  that  the  segment  has  fixed  state  and  time  horizon 
boundaries  and  the  control  and  state  trajectories  will  lie  on  the  original  trajectory.  In  the 
problem  formulation  in  DIDO©  this  poses  a  difficulty  since  to  solve  an  optimal  control 
problem,  one  needs  to  detennine  the  cost  to  be  minimized.  If  both  mass  (part  of  the  state) 
and  time  are  fixed  boundary  conditions,  then  there  is  no  cost  to  minimize.  Therefore,  the 
problem  was  formulated  to  minimize  the  final  fuel  cost,  fixing  the  initial  boundary 
condition  at  some  x0  =  [x°,  y°  ,z° ,vx° ,v°  ,v.°  ,m° and  the  final  boundary  condition  at 

some  xf  =[x/,y/,z/,v/,vJ,/,v/]  such  that  the  initial  and  final  states  are  relatively 

close  to  the  initiation  and  termination  of  the  maneuver.  This  problem  is  solved  with  the 
knowledge  that  an  optimal  solution  can  be  found  on  the  segment,  but  that  the  cost  can  be 
no  less  than  the  segment  on  the  original  optimal  path. 

A  total  of  six  experiments  were  run,  using  three  different  sets  of  boundary 
conditions  each  associated  with  a  different  time  horizon  around  the  original  singular  arc 
(TH1,  TH2,  and  TH3).  For  each  specified  time  horizon,  DIDO©  was  supplied  with  two 
initial  guesses  to  choose  from:  (1)  the  original,  coarse  solution  bound  by  the  time 
horizon,  and  (2)  a  two-point  guess  using  only  the  initial  and  final  nodes  of  the  original 
primal  solution  bounded  by  the  time  horizon.  Figure  51  shows  the  thrust  profile  for  the 
singular  arc  in  the  original  main  engine  solution  and  the  three  time  horizons  examined. 
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Figure  5 1 .  Three  Time  Horizons  Around  Original  Main  Engine  Singular  Arc 


The  spacecraft  main  engine  characteristics  (Scarritt  et  ah,  2009)  and  fixed  initial 
and  terminal  conditions  are  as  follows: 

•  Main  Engine  Thrust:  33,631.6621  N 

•  Main  Engine  Isp:  326  sec 

•  Initial  and  terminal  boundary  conditions  for  TH1 


x0 

-3,130.762702562802  km 

To 

-5,554.427470301326  km 

zo 

23,156.51130870926  km 

TcO 

= 

-0.073217640485400  km  / sec 

v,0 

0.189457232407054  km  / sec 

E-o 

-0.246837168808192  km/ sec 

_mo  _ 

16,799.81436532424  kg 

xf 

-3,069.163608399203  km 

y f 

-5,150.061768188334  km 

zf 

22,667.78130088386  km 

= 

0.173269593413525  km/ sec 

Vyf 

0.288498099224656  km/ sec 

V 

-0.324760205524167  km  /sec 

m. 

Free 

L  J 
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Initial  and  terminal  boundary  conditions  for  TH2 


x0 

-3,098.602243204757  km 

Vo 

-5,637.287708143786  km 

zo 

23,263.91379554111  km 

V,0 

= 

-0.073590413655685  km  /  sec 

v,° 

0.188562680976605  km  / sec 

V-0 

-0.243173237675902  km  /  sec 

_m0  \ 

16,799.81436532424  kg 

1 

Xf 

-2,995.8032622945048  km 

yf 

-5,026.665289857207  km 

zf 

22,528.09830807556  km 

v*/ 

= 

0.169676878511166  km/ sec 

v>/ 

V=f 

0.287787468988926  km  /  sec 

-0.327515426089892  km/ sec 

Free 

mf 

L  J 

•  Initial  and  terminal  boundary  conditions  for  TH3 


-3,066.151333338768  km 

Vo 

-5,720.134711415117  km 

zo 

23,370.19755922657  km 

v,o 

= 

-0.073812293999072  km  /sec 

V,0 

0.187756011321829  km  /sec 

V-0 

-0.239609189112516  km /sec 

_mo  J 

16,798.21977554951  kg 

xf 

-2,923.002666940919  km 

V/ 

-4,903.646847991885  k/u 

Z/ 

22,387.55942551430  km 

v*/ 

= 

0.171736936444797  km/ sec 

Vyf 

0.289275403273065  km  /sec 

V 

-0.331764528528468  km  /sec 

Free 

»2/ 

L  J 

Theoretically,  the  initial  mass  for  each  of  time  horizon  should  be  the  same,  since 
from  the  original  solution  the  thruster  is  not  firing,  thus  not  burning  fuel  at  the  first  three 
nodes  of  TH3.  The  differences  come  from  the  sensitivity  of  the  scaled  mass  values.  The 
scaled  initial  mass  value  for  TH3  starts  to  differ  at  the  fourth  significant  digit  and  results 
in  a  difference  of  about  1.6  kg. 

•  TH1  and  TH2  Initial  Mass  =  0.825953636218675 

•  TH3  Initial  Mass  =  0.825875239089155 
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A. 


PROBLEM  FORMULATION 


The  problem  formulations  for  each  time  horizon  are  identical,  with  only  the  initial 
and  terminal  boundary  conditions  changed  for  each  case.  For  each  case,  two  different 
initial  guesses  were  provided,  which  in  each  case  influenced  the  final  solution. 

Minimum  Time  Formulation 
Where  x  =  [x,  y,  z,  vx,  vv,  vz,  m]  and  u  =  [T,a,/3], 


Minimize 
Subject  to 
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B.  TH1  RESULTS 

Solving  this  problem  using  the  TH1  boundary  conditions  and  supplying  two 
different  initial  guesses  to  DIDO©  results  in  two  different  solutions,  each  satisfying  the 
necessary  conditions  for  optimality.  The  solutions  are  each  solved  up  to  80  LGL  points 
each.  Comparing  the  two  cases,  the  control  and  state  trajectories  differ  between  the  fixed 
boundary  conditions,  however  both  use  the  same  amount  of  fuel.  Figures  52,  53,  and  54 
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show  the  resulting  control  trajectories  compared  to  the  initial  guess  trajectory  provided. 
Note  that  the  coarse  guess  is  the  same  trajectory  given  by  the  original  problem  using  the 
main  engine.  Also,  the  differences  in  the  fuel  consumed  are  negligible  as  shown  in 
Table  6.  The  solutions  for  each  of  the  time  horizons  are  compared  against  the  original 
coarse  solution  for  the  corresponding  time  horizon.  The  fuel  consumed  for  each  region 
has  differences  due  to  numerical  variations  of  the  scaled  mass  in  the  original  coarse 
solution,  therefore  the  fuel  consumed  in  each  of  the  regions  appear  to  have  different 
values.  In  scaled  units  the  differences  are  small  between  the  time  horizons. 


Fuel  Consumed  (kg) 

Variation 

Original  Coarse  Solution 

1,360.2 

— 

Using  Coarse  Guess 

1,361.0 

0.0006  % 

Using  2-point  Guess 

1,361.2 

0.0007  % 

Table  6.  TH1  Fuel  Consumption 


Figure  52.  Thrust  Trajectories  for  TH1 
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P  (rad) 


Figure  53.  Azimuth  Angle  Trajectories  for  TH1 


Figure  54.  Elevation  Angle  Trajectories  for  TH1 
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Of  note  is  the  aliasing  that  occurs  when  the  initial  coarse  solution  is  used  as  an 
initial  guess  (Ross,  I.M.,  Gong,  Q.,  and  Sehkavat,  P.,  2007).  The  higher  fidelity  solution 
follows  closely  with  the  original  solution.  Similarly,  the  state  trajectories  are  also  closely 
aliased  when  given  the  original  solution.  Using  the  two-point  guess,  the  state  and  control 
trajectories  deviate  from  the  original  solution;  however,  the  final  minimized  cost  is  the 
same.  Therefore,  both  of  these  solutions  appear  feasible  and  as  described  below  they 
both  meet  the  necessary  conditions  for  optimality. 

Figures  55,  56,  and  57  show  that  the  Hamiltonian  is  constant,  the  end  point 
transversality  condition  is  met,  and  the  primer  angle  is  equal  to  180  degrees  and  thus,  the 
Hamiltonian  Minimization  Conditions  are  met.  Additionally,  both  solutions  meet  the 
necessary  conditions  for  optimality  for  singular  arcs  and  the  Legendre-Clebesch 
condition  as  shown  in  Figures  58-63. 


Figure  55.  Hamiltonian  Evolution  (TH1) 
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Figure  60.  Fourth  Time  Derivative  of  Switching  Function  for  Coarse  Guess  Solution 


Figure  6 1 .  Fourth  Time  Derivative  of  Switching  Function  for  2-Point  Guess  Solution 

88 


Figure  62.  Generalized  Legendre-Clebesch  Condition  for  Coarse  Guess  Solution 


Figure  63.  Generalized  Legendre-Clebesch  Condition  for  2-Point  Guess  Solution 
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Examination  of  the  state  vector  trajectories  shows  that  each  of  the  two  cases  for 
the  TH1  bounds  start  and  complete  the  maneuver  at  the  specified  state  boundary 
conditions.  While  the  final  mass  was  not  posed  as  a  constraint  on  the  problem,  the 
resulting  final  mass  is  the  same  as  in  the  original  main  engine  problem  within  very  tight 
tolerances.  In  Figures  64  and  65,  the  position  and  velocity  state  vector  trajectories  for  the 
TH1  singular  arc  problem  are  shown  and  are  described  below.  The  position,  or 
displacement,  state  trajectories  differ  slightly,  however  the  velocity  state  trajectories  are 
more  varied  due  to  the  immediate  relationship  between  thrust  and  change  in  velocity. 


r  =  yj x 2  +  y2  +  z2 
V  =  ^vv2  +  v,2  +  V  2 
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Figure  65.  Velocity  State  (v  )  Trajectories  For  TH1 


The  overall  change  in  mass  is  the  same  for  the  two  feasible  solutions  as  shown  in 
Figure  66.  The  mass  state  trajectory  is  clearly  related  to  the  shorter  duration,  but  higher 
magnitude  of  thrust  in  the  2-point  guess  solution  as  compared  to  the  solution  from  the 
coarse  guess. 
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Figure  66.  Fuel  Consumption  for  TH1 


C.  TH2  RESULTS 

The  same  analysis  conducted  for  the  TH1  solutions  was  performed  for  TH2  and 
TFI3.  The  Flamiltonian  Minimization  Conditions  for  optimality  were  met,  as  well  as  the 
additional  conditions  for  singular  arcs.  Each  case  resulted  in  different  feasible  solutions 
for  the  controls,  with  the  state  variable  trajectories  satisfying  the  fixed  boundary 
conditions.  And  while  the  state  trajectories  vary  slightly  between  the  cases,  the  fuel 
consumption  is  the  same  between  the  original  coarse  solution  and  the  optimal  solutions. 
These  results  further  indicate  the  range  of  possibilities  for  conducting  the  singular  arc 
maneuver. 

In  TH1,  the  maximum  thrust  exceeded  the  maximum  thrust  of  the  auxiliary 
engines;  however  the  resulting  control  trajectories  for  TFI2  and  TFI3  achieved  maximum 
thrust  less  than  the  auxiliary  engines  (4.4  kN).  This  gives  some  indication  that  the 
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auxiliary  engines  can  indeed  be  used  for  the  singular  burn  maneuver;  however  the  fuel 
consumption  will  necessarily  be  greater  due  to  the  lower  rocket  specific  impulse  as 
compared  to  the  main  engines. 

Figures  67,  68,  and  69  show  the  control  trajectories  for  TH2.  The  maximum 
thrust  for  both  solutions  is  about  the  same,  around  4  kN.  The  control  angles  during  the 
bum  are  identical  to  the  angles  arrived  at  in  the  TH1  solutions.  They  are  also  found  to  be 
equal  in  the  TH3  solutions  as  well.  Clearly,  a  rocket  engine  with  maximum  thrust  of  4.4 
kN  would  be  able  to  execute  these  maneuvers.  The  difference  will  be  in  fuel 
consumption  alone. 


Figure  67.  Thrust  Trajectories  for  TH2 
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P  (rad) 


Figure  68.  Azimuth  Angle  Trajectories  for  TH2 


Figure  69.  Elevation  Angle  Trajectories  for  TH2 
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Displacement  State  Trajectory  (scaled) 


Neither  of  the  state  trajectories  have  the  same  aliasing  that  was  seen  in  TH1, 
however  the  two  solutions,  using  different  initial  starting  guesses,  are  very  similar  as  seen 
by  the  control  trajectories  above  and  the  state  trajectories  shown  in  Figures  70,  71  and  72. 


Figure  70.  Displacement  State  (r)  Trajectories  for  TH2 
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Figure  7 1 .  Velocity  State  ( v  )  Trajectories  for  TFI2 


Figure  72. 


Fuel  Consumption  for  TH2 
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The  fuel  consumption  for  TH2  is  listed  in  Table  7  and  the  differences  are  again 
negligible  as  compared  to  the  original  coarse  solution  within  the  specific  time  horizon. 


Fuel  Consumed  (kg) 

Variation 

Original  Coarse  Solution 

1,339.7 

— 

Using  Coarse  Guess 

1,339.1 

0.0007  % 

Using  2-point  Guess 

1,339.2 

0.0004  % 

Table  7.  TH2  Fuel  Consumption 


D.  TH3  RESULTS 

By  expanding  the  time  horizon  to  TH3,  the  control  trajectories  begin  to  deviate 
from  the  original  solution.  These  solutions  also  meet  the  optimality  conditions  for 
singular  arcs  and  are  also  feasible  solutions.  In  both  solutions  arrived  at  from  different 
initial  guesses,  the  maximum  thrust  is  again  less  than  4.4  kN  indicating  that  this  particular 
control  profile  could  be  used  with  the  auxiliary  engines.  The  solution  from  the  2-point 
guess  is  interesting  in  that  the  thrust  profile  is  starting  to  appear  as  a  bang-bang 
maneuver,  however  in  this  situation  remains  a  singular  arc  solution.  In  general,  these 
solutions  are  not  as  “clean”  from  a  control  perspective  because  neither  solution  has  a  zero 
to  zero  thrust  trajectory.  In  both  cases  there  is  a  non-zero  initial  thrust  as  seen  in  Figure 
73.  In  Figures  74  and  75  the  control  angles  are  approximately  the  same  for  both  initial 
guesses  even  though  the  coarse  gain  has  a  single  burn  while  the  2-point  guess  solution 
has  a  two  burn  solution. 
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a  (rad)  Thrust  (KN 


Figure  73.  Thrust  Trajectories  for  TH3 


Figure  74.  Azimuth  Angle  Trajectories  for  TH3 
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Figure  75.  Elevation  Angle  Trajectories  for  TH3 


The  state  trajectories  as  shown  in  Figures  76,  77,  and  78  illustrate  a  more  extreme 
example  of  varying  state  trajectories  with  the  minimized  cost  of  fuel  being  equal, 
emphasizing  the  result  that  multiple  local  minimal  solutions  exist  around  the  singular  arc. 
Table  8  shows  the  negligible  difference  in  fuel  from  the  original  coarse  solution  for  the 
given  time  horizon. 
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Displacement  State  Trajectory  (scaled) 


Figure  76.  Displacement  State  Trajectories  (TFI3) 


Figure  77. 


Velocity  Trajectories  (TH3) 
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Figure  78.  Fuel  Consumption  (TH3) 


Fuel  Consumed  (kg) 

Variation 

Original  Coarse  Solution 

1,345.8 

— 

Using  Coarse  Guess 

1,346.3 

0.0007  % 

Using  2-point  Guess 

1,346.7 

0.0004  % 

Table  8.  Fuel  Consumption  for  TH3 


E.  FEASIBILITY  OF  USING  AUXILIARY  ENGINES  FOR  SINGULAR  ARC 
MANEUVER 

In  the  analysis  of  the  singular  arc  using  the  main  engine,  it  was  found  that  for  the 
TFI2  and  TFI3  time  horizons,  it  appeared  feasible  to  use  the  auxiliary  engines  in  place  of 
the  main  engines  to  conduct  the  singular  arc  maneuver.  However,  the  TH1  solutions 
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using  the  main  engine  resulted  in  thrust  profiles  that  exceeded  the  capacity  of  the 
auxiliary  engines.  It  was  expected  that  the  lower  Isp  would  result  in  a  higher  fuel 
consumption  as  compared  to  the  main  engine,  given  the  same  magnitude  of  thrust.  The 
same  six  simulations  as  examined  in  the  previous  section  were  perfonned  using  the  same 
optimal  control  problem  formulation,  but  using  the  thrust  and  Isp  parameters  for  the 
auxiliary  engines. 

As  expected,  the  control  profiles  around  TH2  and  TH3  using  the  auxiliary  engines 
were  identical  to  those  found  using  the  main  engines.  The  TH1  solution  for  the  auxiliary 
engines  however,  was  constrained  in  maximum  thrust  and  resulted  in  a  different,  yet 
feasible  control  structure.  Figure  79  shows  the  different  thrust  profiles  attained  in  the 
TH1  time  horizon  using  the  same  initial  primal  guess  from  the  main  engine  singular  arc. 
In  this  case,  the  two  different  initial  guesses  resulted  in  very  similar  control  profiles. 


Figure  79.  Auxiliary  Engine  Thrust  Profile  (TH1) 
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In  all  three  time  horizons,  using  the  auxiliary  engines  resulted  in  an  increase  in 
fuel  consumption  as  shown  in  Table  9.  The  results  show  that  for  a  relatively  small 
increase  of  fuel,  the  singular  arc  maneuver  can  be  conducted  using  the  auxiliary  engines. 
In  order  to  solve  for  the  total  trajectory  fuel  consumption  where  the  main  engines  are 
used  for  the  first  and  third  TEI  burn  and  the  auxiliary  engine  is  used  for  the  singular  arc 
maneuver,  another  step  using  the  minimum  fuel  problem  formulation  needs  to  be 
conducted  between  the  boundary  conditions  of  the  singular  arc  time  horizon  and  the 
original  boundary  conditions  provided  to  solve  the  full  moon-to-earth  trajectory  problem. 


Main  Engine 

(33.1  kN) 

Aux  Engine 

(4.4  kN) 

Singular  Arc  on  Full  Trajectory  (Yan  et  ah,  2010) 

1,333.6  kg 

— 

TH1 

1,361.0  kg 

1,432.6  kg 

TH2 

1,339.1  kg 

1,409.6  kg 

TH3 

1,346.3  kg 

1,417.1  kg 

Table  9.  Fuel  Consumption  Comparison  Between  Main  Engine  and  Auxiliary  Engines 

Over  Singular  Arc 


F.  CONCLUSIONS 

The  reason  that  multiple  feasible  and  optimal  solutions  appear  is  because  of  the 
existence  of  the  singular  arc.  Recalling  that  the  switching  function  equals  zero  for  a 
singular  arc,  it  is  clear  that  the  resulting  magnitude  of  thrust  can  lie  between  the  minimum 
and  maximum  thrust.  If  the  initial  and  final  boundary  conditions  are  established  at 
varying  positions  in  time  before  and  after  the  singular  arc,  it  is  possible  to  generate 
different  control  solutions  and  thus,  different  state  trajectories  that  are  also  optimal 
solutions.  This  gives  rise  to  the  possibility  of  multiple  feasible  solutions  between  two 
given  boundary  conditions  with  a  singular  arc  between  them.  From  an  engineering  point 
of  view,  this  offers  greater  flexibility  in  designing  a  solution  that  closer  resembles  the 
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imposed  physical  constraints  on  the  available  controllers.  The  difficulty  though,  is  in 
finding  a  fuel  optimal  solution  that  meets  the  given  constraints  and  requirements  from  a 
possible  myriad  of  solutions.  For  example,  while  some  modern  rocket  engines  may  have 
the  ability  to  throttle,  many  typical  rockets  are  forced  to  apply  a  bang-bang  type  control. 
In  seeking  for  a  bang-bang  solution  to  perform  the  same  maneuver  as  the  singular  arc, 
additional  constraints  will  have  to  be  imposed,  thus  changing  the  problem  formulation.  It 
is  possible  to  use  pseudo-spectral  tools  such  as  DIDO©  to  find  a  range  of  optimal  bang- 
bang  solutions,  however  it  is  likely  that  the  fuel  cost  will  be  greater. 

The  singular  arc  maneuver  results  in  a  slight  change  in  the  orbits  eccentricity 
(semi-major  axis);  however  it  appears  that  the  maneuver  is  used  primarily  to  change  the 
orbit  plane  (inclination).  Further  study  of  orbit  maneuvers  may  show  that  singular  arcs 
will  appear  in  fuel  optimal  orbit  transfers  whenever  an  inclination  change  occurs, 
assuming  that  the  rocket  thrust  is  allowed  the  full  range  from  null  to  maximum.  To 
further  examine  the  nature  of  fuel-optimal  plane  change  maneuvers,  it  would  be  prudent 
to  examine  pure  plane  change  maneuvers  as  well  as  combined  plane  change  and  apogee 
raising  (or  lowering)  maneuvers  using  similar  n-body  dynamics  equations  as  used  here. 
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APPENDIX  A.  SINGULAR  ARC  NECESSARY  CONDITIONS 


From  Park  et  al.,  2010: 


Third  derivative 
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Fourth  derivative 


=  -\  ‘  [5,]+  |  _  —[5-3] 
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[53],  where[Si]  =  YjSi 
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J4,4 
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Legendre-Clebesch  Condition  (convexity) 


d  d4S 
dT  dt4 


>  0 


APPENDIX  B.  END-POINT  TRAN S VERS ALIT Y  CONDITION 


The  transversality  condition  states: 

A( tf)  = - ,  where  E{v_,x[tf^  =  istx((/-))  +  ifre  (*(*/)) 

dXf 

The  variable,  v ,  represents  the  Lagrange  multiplier  vector  associated  with  the  end- 
point  state  condition, e^x [t f )) .  The  end-point  cost  isisfx^n.  The  given  end-point 

boundary  condition  is  defined  asxy  =  [x7,  y  f ,  zf ,  v/ ,  v/ ,  v/] .  As  can  be  seen  in  each 
of  the  problems,  the  only  additional  information  yielded  is  Xm  (  tf  j . 

Minimum  Fuel  Problem  End-point  Transversality  Condition 
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Minimum  Time  Problem  End-point  Transversality  Condition 
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